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ELECTROMAGNETIC  ENERGY  DEPOSITION  IN  A  CONCENTRIC 
SPHERICAL  MODEL  OF  THE  HUMAN  OR  ANIMAL  HEAD 

INTRODUCTION 

The  head  is  modeled  by  several  homogeneous  regions  of  tissue  bound¬ 
ed  by  one  or  two  members  of  a  family  of  concentric  spheres.  These  tis¬ 
sues  include  brain,  cerebrospinal  fluid  (CSF),  dura,  bone,  subcutaneous 
fat,  and  skin  tissue.  We  assume  that  this  complex  of  biological  material 
is  exposed  to  nonionizing  electromagnetic  radiation  taking  the  form  of  a 
time-harmonic  plane  wave  of  peak  amplitude,  Eq.  The  time  variation 
factor,  exp(-iwt),  has  been  suppressed  in  most  of  the  discussion. 

Wave  propagation  is  in  the  positive  z-direction,  and  the  electric  field, 
E,  is  linearly  polarized  in  the  x-direction  (Fig.  1).  A  rectangular- 
spherical  coordinate  system  with  origin  at  the  center  of  an  inner  core 
sphere  is  used.  Also,  the  medium  surrounding  the  concentric  spherical 
model  is  taken  as  free  space  (or  vacuum).  Thus  our  embedding  medium  is 
a  nonconductor,  and  both  the  surrounding  medium  and  the  model  are  non¬ 
magnetic.  Each  region  (p  =1,...,N-1)  into  which  the  model  is  parti¬ 
tioned  is  homogeneous,  isotropic,  and  possesses  a  unique  dielectric 
constant  and  conductivity.  All  magnetic  permeabilities  are  considered 
to  have  the  value  unity.  The  value  "p  =  N"  is  reserved  for  reference 
to  the  containing  medium. 

The  need  for  a  multilayer  model  and  the  inadequacy  of  (1)  ignoring 
the  relatively  thin  outer  structures,  or  (2)  carrying  out  a  volume 
average  of  the  electrical  properties  of  the  regions  can  be  seen  by 
looking  at  Figure  2  for  the  case  of  the  rhesus  monkey.  There  graphi¬ 
cally  displayed,  for  comparison  purposes,  are  three  superimposed  dis¬ 
tributions  of  absorbed-power  density  along  the  z-axis.  The  monkey-head 
models  consist  of  (1)  pure  brain  tissue,  (2)  tissue  with  average  volume 
of  electrical  properties  of  the  structural  components  in  Figure  1,  and 
(3)  unique  tissues  represented  in  Figure  1.  The  predicted  distributions 
are  based  on  1-volt-per-meter  intensity  incident  wave  at  3  GHz. 
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FTgure  1.  Electromagnetic  plane  wave  impinging  on  a  head 
model  composed  of  an  inner  core  sphere  and  five 
spherical  shells. 

Dimensions  of  structural  media  and  electrical  parameter  values  were 
extracted  from  a  table  that  was  produced  by  Shapiro  et  al .  (13).  Table 
1  presents  such  information.  Contour  plots--Figures  3  (<j>  =  0) ,  4  (<|>  = 
tt/2),  and  5--are  likewise  based  on  information  offered  by  Table  1. 

The  linear  plot.  Figure  6,  and  the  contour  plots--Figures  7  (<)>  =  0), 
8  (<j>  =  tt/2),  and  9--are  founded  on  the  entries  of  Table  2,  an  extraction 
from  a  paper  by  Weil  (15).  Incident  plane-wave  characteristics  are 
1-volt-per-meter  intensity  and  1-GHz  frequency.  Other  parameter  values 
pertinent  to  the  computations  for  graphical  construction  are  given  in 
Table  2. 


POWER  DENSITY  ( WRT TS  PER  CUBIC  METER) 


TABLE  1.  RHESUS-MONKEY- HEAD  DATA 


Region  Tissue  Thickness  Radius  of  Relative  Conductivity,3 

(p)  modeled  (cm)  surface  dielectric  o 

boundary,  cortsta„t>o  (inhJ/m) 

.  P.  £n 


1 

Brain 

sphere 

2.68 

42.0 

2.0 

t. 

2 

CSF 

0.20 

2.88 

77.0 

1.9 

3 

Dura 

0.05 

2.93 

45.0 

2.5 

4 

Bone 

0.20 

3.13 

5.0 

0.2 

i 

5 

Fat 

0.07 

3.20 

5.0 

0.2 

6 

Skin 

0.10 

3.30 

45.0 

2.5 

% 

* 

aAt  T  =  37 

°C  and  f 

=  3  GHz. 

: 

TABLE  2. 

IDEALIZED  HUMAN-HEAD  DATA 

Region 

Tissue 

Thickness 

Radius  of 

Relative 

Conductivity,3 

(P) 

modeled 

(cm) 

surface 

dielectric 

boundary, 

constant,3 

P 

(mho/m) 

rp 

E. 

(cm) 

P 

1 

Brain 

sphere 

9.10 

60.00 

0.90 

2 

CSF 

0.20 

9.30 

76.00 

1.70 

3 

Dura 

0.05 

9.35 

45.00 

1.00 

4 

Bone 

0.40 

9.75 

8.50 

0.11 

5 

Fat 

0.15 

9.90 

5.50 

0.08 

6 

Skin 

0.10 

10.00 

45.00 

1.00 

aAt  f  =  1  GHz. 
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POWER  DENSITY 
(W/M3) 

H  0.295 
Y  0.265 
Z  0.235 
X  0.205 
*  0.175 
O  0.145 
X  0.115 
+  0.084 
a  0.055 
O  0.025 


Figure  3.  Distribution  of  power  density  in  a  3.3-cm-radius 
sphere,  a  model  of  the  rhesus-monkey  head,  at 
3  GHz.  (E-plane) 


7 


Figure  4.  Distribution  of  power  density  in  a 
3. 3-cm- radius  sphere,  a  model  of  the 
rhesus-monkey  head,  at  3  GHz.  (H-plane) 
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POWER  DENSITY 
(W/M3) 

K  0.0595 
Y  0.0535 
Z  0.0475 
X  0.0415 
*  0.0355 
O  0.0295 
X  0 .0235 
4  0.0175 
A  0 .01  15 
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Nor.uniform  distribution  of  the  absorbed  energy  inside  of  spheres 
gives  rise  to  the  internal  appearance  of  "hot  spots."  Kritikos  and 
Schwan  (5)  showed  that  "hot  spots"  appear  in  lossy,  homogeneous,  brain¬ 
like  material  spheres  for  the  range  of  values  of  the  radius  and  frequen¬ 
cy:  0.1<R<8  cm  and  300<f<1200  MHz.  Shapiro  et  al .  (13)  and  Weil  (15) 

have  shown  the  existence  of  hot  spots  inside  multilayered  spheres. 

Each  has  taken  into  account  the  importance  of  the  frequency  dependence 
of  the  sphere  electrical  properties.  What  precise  conditions  will  pre¬ 
cipitate  the  phenomenon  are  still  not  well  known.  We  do  know  that  the 
radius  cf  the  sphere  and  frequency,  among  others,  do  play  a  significant 
role.  Occurrence  takes  place  at  the  front  surface  or  inside  the  head. 

It  is  a  phenomenon  that  is  importantly  connected  to  small  animals  and 
infants . 

The  concentric  spherical  model  represents  one  step  forward  in 
approximating  reality  as  compared  to  the  single,  homogeneous  sphere. 

Even  so,  the  shortcomings  of  this  model  are  to  be  found  in  (1)  shape, 

(2)  electrical  properties,  (3)  thicknesses  of  tissue  media,  (4)  assump¬ 
tion  of  tissue  media  being  homogeneous  and  isotropic,  and  (5)  inoperative 
conduction,  convection,  and  radiation-heat-transfer  mechanisms.  The 
present  computer  program  will  be  updated  in  the  near  future  by  an 
attempt  to  incorporate  the  mechanism  of  blood  flow  (along  with  other 
features).  This  may  result  in  an  appreciable  reduction  in  the  temper¬ 
ature  rise  now  calculated. 

The  knowledge  to  be  gleaned  from  this  current  research  is  directly 
related  to  the  research  effort  of  the  Radiation  Sciences  Division  at  the 
USAF  School  of  Aerospace  Medicine.  Briefly,  here  studies  are  being  con 
ducted  to  (1)  determine  the  radiofrequency  radiation-induced  effects  in 
biological  specimens,  (2)  seek  out  possible  hazards  to  personnel  in  a 
radiofrequency  environment,  (3)  accurately  measure  and  determine  the 
distribution  of  energy  In  the  whole  biological  body  or  just  in  a  partic¬ 
ular  organ,  (4)  extrapolate  response  to  radiation  from  the  test  animal 


to  man  in  a  meaningful  manner,  and  (5)  contribute  in  the  design  of 
realistic  safety  standards  with  a  solid  basis. 

The  division  of  this  paper  entitled  "Mathematical  Description" 
consists  of  four  sections.  Since  spherical  harmonics  (Stratton  [14], 
pp.  399-423)  are  used  to  expand  the  incident,  induced,  and  scattered 
fields,  we  include  In  the  section  "Mathematical  Preliminaries"  details 
of  the  exact  evaluation  of  inner  products  entering  into  the  computa¬ 
tion  of  expansion  coefficients.  In  the  section  "Expansion  of  Induced 
Fields  in  Terms  of  Vector  Wave  Functions,"  the  expansions  are  used  to 
solve  Maxwell's  equations  (Stratton  [14],  p.  26),  subject  to  the  con¬ 
dition  that  the  tangential  components  of  electric  field  t  and  magnetic 
field  fi  are  continuous  across  the  spheres  delimiting  different  regions 
of  the  head  model.  The  section  "Determination  of  Total  Absorbed  Power" 
considers  the  integrals  that  appear  in  Poynting's  theorem.  Such  Inte¬ 
grals  are  evaluated  in  closed  form,  thereby  yielding  a  formula  for  the 
total  power  absorbed.  The  section  "Summary  of  Key  Equations  and 
Formulas"  contains  a  detailed  summary  of  the  formulas  Implemented  by 
the  computer  program.  Concentric  Spherical  Model  (CSM),  for  automati¬ 
cally  calculating  the  radiofrequency  electromagnetic  energy  deposited 
in  a  concentric  spherical  model  of  the  human  or  animal  head. 

The  succeeding  division,  entitled  "Program  Description,"  contains 
pertinent  information  about  the  computer  program  that  is  beneficial 
to  the  user.  The  appendixes  consist  of  a  sample  problem  with  computer 
results  and  a  source  listing  of  the  program  CSM. 

To  benefit  users  of  this  report,  program  CSM  is  described  in 
sufficient  depth  to  permit  job  setups  and  Implementation  on  any  modern 
computer.  The  mathematical  theory  and  formulas  basic  in  accomplishing 
the  computations  are  discussed  in  an  extensive  manner.  Discussion  of 
this  users-oriented  computer  program  covers  such  details  as  structure 
and  sequence  of  control  parameter  and  data  cards,  output  formats,  and 
subroutine  and  function  subprograms.  Sample  printouts  and  plots  (linear 
and  contour)  of  computer  results  and  a  listing  of  the  FORTRAN  IV  source 
program  are  included. 
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MATHEMATICAL  DESCRIPTION 


Mathematical  Preliminaries 


This  part  of  the  mathematical  discussion  introduces  an  inner  prod¬ 
uct  on  doubly  oeriodic  vector  valued  functions,  presents  a  study  of 
some  of  the  properties  of  Legendre  polynomials  used  in  evaluating  inner 
products,  introduces  vector  wave  functions,  and  verifies  some  of  their 
properties. 

Definition  1.  Let  S  =  {(8,<f>):  0  <_  9  <  ir  and  0  <J>  <  tt}.  Then 

let  H ( S , C )  denote  the  continuous  complex  valued  functions  A  defined 
on  S  that  satisfy  the  inequality 


•2ir  ('tt  ^ 

[  |A(e,4>)|  sinede]d4>  <  “  . 

Jo  'o 


(1) 


For  any  functions ,  A  and  8  in  H(S,C)  define  the  inner  product 
<,>  by  the  rule 


<  A,B  > 


■2  IT  j'TT  - 

[  A(0,<ji)B(e,<j>)sinede]d()>  . 

o  o 


(2) 


Rroposition  1.  The  space  H(S,C)  with  the  inner  product  <»>  is 
a  pre-Hilbert  space. 

This  follows  immediately  from  the  definition. 

Now  we  need  some  properties  of  the  associated  Legendre  functions. 

Definition  2.  For  alt  nonnepa! :  <’. •  ir.tearvs,  m  and  n  define 


2  m/2 

p>)  •  ik»_>  D"+V-i)n . 


(3) 
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Proposition  2.  If  m+n  is  even  (odd),  then  Dn+m(x2-l)n  is  a 
linear  combination  of  even  (odd)  pou'ers  of  x. 

Proof.  Observe  that 

(x2-l)"=  l  (")x2k(-U"-k  .  (4) 

k=o 

Since  an  even  (odd)  number  of  derivatives  of  an  even  power  of  x  is  an 
even  (odd)  power  of  x,  the  proposition  is  true. 

Corollary  1.  If  n-Hn  is  even  (odd),  then  p^(x)  is  an  even  (odd) 
function  of  x . 

Proof  of  corollary  1.  Since  (l-x2)m^2  is  an  even  function  of  x, 
the  corollary  follows  immediately  from  Proposition  2. 

Proposition  3.  For  all  nonnegative  integers  m  and  n 


I  = 


Ocos0)^(Ocose))sin0d0 


0  . 


Proof.  Since 

^d?)(Pn(cos0)2)  =  Ocos0)^)(Ocose))  ’ 


(5) 


(6) 


an  integration  by  parts  implies  that 


I  =  ^[P^cose^sinejg 


P^(cos9)2cosed0]  . 


(7) 
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Substituting  x  =  cose  and  using  the  fact  that 


de  =  -  (1/ 


)dx 


(8) 


it  follows  that 


I  « 


% 


A 

P%)2(x/  /l-?  )dx  , 
"i  n 


(9) 


which  in  view  of  Corollary  1,  implies  that  I  =  0. 

Proposition  4.  For  c.11  nonnegative  integers  m  and  n,  we  have 

that 


P™(cose)P™(cose)sin9de 

o 


fO.r^n 

2(n+m)!  ,  r=n. 

L(2n+l)(n-m)» 


(10) 


Proof.  The  proof  is  carried  out  completely  in  Whittaker  and  Watson 
(16,  pp.  324-325). 

Proposition  5.  For  all  nonnegative  integers  m,  n,  and  r,  we 

have 


^(P^cose))^-(p”(c°s6))sine<le  •  A"n>r) 


(--*-)  p'''(cose)P'"(cose)sinede. 
o  sin  e  n 


(11) 


Proof.  Observe  that 


in  n  m/2'1 

A  /,  '  ,n+m  n  ^(1-x  ;  /  n  N  ,n+m  ,  2  ,Nn 

dr(l-x  )  _d  ,  2  =  2 v  '  (-2x)  _d_  (x  -1) 

dx  2"n!  dxn+m  2nn!  dxntm 


+  (1-0 


m/2 


2nn ! 


.  n+m+1  0 

—  (x2-l)n 

dxn+m+1  ( 


(12) 


Let  (d/de)F(cose)  =  F' (cose)(-sin8).  Thus,  x  =  cose  implies  that 


dx  d _ d_ 

de  dx  “  de  *  (13) 


Hence 


A(n,r)  =  (  (-sin0)pn  (x) (-sine )pJJ  '(x)sinede 
-  ^1(l-x2)Pj,(x)Pj,(x)(-dx) 

=  4  (l-x2)p^(x)p;'(x)dx  .  (14) 

J_1  n  r 
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Integrating  by  parts  in  the  above  integral,  we  find  that 


W) 


[n(n+l) - yy]  P„(x)P"l(x)dx 

-1  l-x  n  r 


1  m2 

-m 


Pyy- +  r(r+l)]P™(x)P™(x)dx 
-1  1-x  "  r 


(15) 


Since  =  A^p)  the  above  relation  shows  that  if  r  J  n. 


Pp(x)Pp(x)dx  =  0. 


-1 


(16) 


Substituting  back  x  =  cose,  we  find  that 


(n,r)  5(n,r)  2n+l  fn-m) ! 


2  (n+m)! 


n(n+l) 


fir  2 

m  vnm 


(  2  ) Pp(cos8 )P' (cose)sinede 

o  sin^e  n 


(17) 
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Proposition  6.  For  all  nonnegative  integers  m,  n,  and  r,  we 

have 


f  [^-(P^cose))  ^-(P^(cose))  +  (m2/sin20)pjj(cose)pj(cose)]sinede 

'  n 


=  $ 


G 


(n,r)v2n+l 


r>  "<"n>  •  <18> 


Proof  of  Proposition  6.  From  Proposition  5  we  deduce  equation  18. 

■> 

Definition  3.  Let  i,  j,  and  k  denote  the  unit  vectors  in  the 
Cartesian  coordinate  system.  Define 

-*•  -+ 

er  =  sin0cos<i>i  +  sinosin^j  +  cosok,  (19) 


e.  =  cos0cos4>i  +  coses i n<^ j  -  sin0k, 
0 


and 


e,  =  -sin<j>i  +  cos<j>j  . 


(20) 


(21) 
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Def  '.nition  4.  If  S  is  a  surface  in  R  bounded  by  a  simple 

"*■  _  I 

closed-c  vrve  C,  and  A  is  a  C  vector  field  defined  in  a  neighbor¬ 
hood  of  S,  then  curl  (A)  is  a  vector  field  such  that 


-y  f-y-y 

curl(A)-Nda  =(|)A*Tds 
C 


(22) 


where  N  and  T  are,  respectively ,  the  unit  normals  and  the  unit 
tangents  of  S  and  C. 

-4- 

Prcoosition  7.  If  A  is  a  vector  valued  function  of  r ,  6,  and 
$ ,  they. 


cur!  (A)  =  ?4H?[0/89)(sin6A^)  -  (3/3<(>)A0]er 
+  ^[(l/sin9)(3/3<j))Ar-(9/3r)(rA^)]e0 

+  7  [(3/3r)(rA0)  -  ( 3/39 )Ar] .  (23) 


Proof.  This  follows  from  Stokes'  theorem  and  the  fact  that  in 

Cartesian  coordinates  x,y,z,  where  F  =  F  i  +  F  j  +  F  k,  the  curl  of 

x  y  z 

vector  F  is  defined  by 

->  3F  3F  -*■  3F  3F  -*• 

curKF)  =  (jf  -  3^)1  +  (jr  -  3^)3 


3F  3F  -*■ 

♦  -  ^)k. 


3x 


3y 


(24) 
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Proposition  8.  In  spherical  coordinates  if  ip  is  a  function  of 
r,  0,  and  <J>,  then 

=  (l/r2H(3/3r)(r2(3/3rh)] 

+  (l/(r2sine))  [(3/39)  (sine  (3/39  )\|>)) 

+  (l/(r2sin2e))(32/3ci>2)i!/  .  (25) 


Proof.  This  follows  from  the  fact  that  in  Cartesian  coordinates 


A*  =  (32/3X2)+  +  (d2/dy2h  +  (32/3Z2)^ 


(26) 


and  the  coordinate  transforms 

x  =  rsin(e)cosU), 

y  =  rsin(e)sin(ij>) ,  (27) 


and 


z  =  rcos(e). 
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Proposition  9. 


For  any 


function  <j>  of  r,  0,  and  4> , 


A  =  4'rer 


(28) 


implies  that 


=  curl (A) 


=  (l/sine)((3/3<i>H)e9-  ( (a/ae )^)e^  .  (29) 


Proof.  This  follows  by  direct  substitution  of  equation  28  into 
equation  23. 

2  .  . 

Proposition  10.  Suppose  4>  is  a  C  function  satvsfying 


+  k  ijj  =  0, 


(30) 


where  k  is  a  complex  number,  then 


M,  =  curl  Ure  ) 

<P  r 


(31) 


and 


N  =  (l/k)curl (M  ) 


(32) 
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imply  that 


=  [{l/(kr))(3/3r)(r2(3/3r)<iO  +  kn|/]er 
+  (l/(kr) )  (32/3/3e)  (nj))e0 

+  (l/(krsine))(3  /3r3<|>)(r^)e  .  (33) 

In  the  next  section  we  work  out  some  consequences  of  this  propo¬ 
sition  when  the  function  ip  is  the  product  of  a  spherical  Bessel 
function,  a  Legendre  polynomial,  and  a  sine  or  a  cosine. 


Expansion  of  Induced  Fields  in  Terms  of  Vector 
Wave  Functions 


In  determining  the  response  to  a  plane  wave  of  a  union  of  regions 
of  dielectric  material  delimited  by  spheres,  we  use  the  vector  wave 
functions  (of  kp  =  complex  propagation  constant  for  the  p-th  layer  of 
the  dielectric  material): 


M 


(e,j)  _ 
(l.n) 


nWzS(kpr>pJ(cos9)sin*  ee 


Zp(k  r)(d/de)(pJ(cose))cos<j>  e  , 


(34) 
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"ain! =  nWzn(kPr)pi(cose)  cos*  ?8 


-  zp(kpr) (d/de )(P* (cose ))sin<j>  , 


+  ^r(a/3r)(rzj(kpr))(d/de)(pJ(cose))cos4>  e 

•  irF5T!«(3/8r)(rzn(kpr))p5(cos0)s’"*  %  • 

n  1  1 


+  ^~(3/3>")(rzj(kpr))(d/de)(pJ(cose))sin<!. 

+  irTW(3/8r)(rzJ(kpr))pS(cose)cos^  % 

p 

zj(p)  =  Jn(p)  =V  ir/2p  0n+^(p), 


and  Jn+P(p)  and  Y  ^(p)  are  the  Bessel  and  Neuman  functions  of 
order  half-an-odd  integer,  respectively. 

Pvoposion  11.  For  all  nonnegative  integers  m  and  n  and  all- 
integers  .j  and  j'  in  {1,2,3},  we  have 


<  M 


(e.j) 

(l.n) 


J(o  ,r) 

(l.m) 


0, 


(41) 


i5(e.J) 


j  (e.j') 

(l.m) 


A*  2  (n+1) ! 

A6(n,m)  2"n+l  (n-1) ! 


M(o,j)  ufo.j') 
M(l,n)  *  M(l,n)  " 


(42) 


where 


A  =  v  r2  z;j(kr)z;J  (kr)  , 


(43) 


j(o,j) 

'(l.n) 


<  Nft'K  ,  N); x #  > 


J(e.j') 

(l.m) 


=  0  , 


(44) 


|(o,j) 

'(l.n) 


<  N)V’^(  ,  N);-v'  > 


J(o.j') 

'(l.m) 


=  B5(n,m)  2n+T  {n-ljl  n(n+1)+cn 


=  < 


w(e,j) 

N(l,n) 


i(e.j') 

(l,m) 


(45) 
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B  =  it  (^7)2(a/3r)(rzj(kpr))(3/3r)(rzj(kpr)) 


(46) 


and 


C  = 


n  2n+l 


2tl_  (n+l)l.  n2(.n+l-)2zj(k  r)zJ'(k  r) 
(n-1 ) !  ,.2„2  VK pr;znlKp  ; 


kpr 


"  M(l,n)  >  N(l,m)  >  "  0  ’ 


<  J(e,j)  Z(e,y)  _  n 


<  "KS  •  >  ■ « • 


(47) 

(48) 

(49) 


(50) 


Proof.  This  follows  from  the  definitions  and  the  facts  that 


[(DP* (cose)) (DP* (cose))  +  — P*(cose)P*(cose)]sinede 


sin  e 


=  6 


(n,m)  2n+l 


n(n*i) 


(51) 


and 


[P*(cos8)D(P*(cose))  +  P*(cose)D(P*(cose))]de  =  0, 


(52) 


where  D  =  d/de  . 
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Now  we  want  to  develop  formulas  relating  the  fields.  Let  us  write 
the  fields  for  the  p-th  region  as 


E  =c  y  *1  2-f+l  f  m(°»1  )  ih  M(e,l) 

ep  "  Eo  pry  [a(f,P)M(i,£)  •  1b(£,p)N(i,£) 


+  a,  m(o>3)  ia  m(s»3)i 

au,pr(i,D  1b(^p)n(i,£)1 


=  __Ef  y  2m.  Th  ^e*1)  +  1(0,1) 

y0co  Eo  ^  irmj  b(£,p)M(l ,£)  1a(£,p)N(l,£) 


e(£>P)”(i;i) +  ia(4p)N(°i:i)] 


■+/ .  .  \  •>/ .  .  \ 

in  terms  of  the  spherical  vector  functions  and  N(i’^) 

Stratton  (14,  p.  564)  for  function  definitions]  and  the  complex  propaga¬ 
tion  constant  kp.  The  tangential  components  of  the  fields  are 


(Ep>«-£„  2*1 


P^(cose) 


p  9’"°  1=1  K7+TT  [aU,p)  — iW”{cos<>)j2(kprp) 

’  ib(£,p)(1/kprp^(9/3r^rJ^kpr)^^d/d0^p^cose)^  cosi’ 

r=rP 

+  _  PJcose)  , 

(tp>  sine  (cos»)h^(kpr!>) 

“  i6(£,p)(1/kprpH(3/9',')(rh^(kpr)))((d/d0)pJ(cose))cos^] ,  (55) 
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(Ep)<j>  "  Eo  ^  ~  1(1+1)  ^a(£,p)('J^kprp)))((d/d0)(P^(cose)))sin<f. 

,  ,  Pp(cose) 

-  ib(£,p)(-1/kprp)((3/ar)(  V(kpr))j(~  ~sinl  )sin* 


P  P  'J.=r  sine 
P 


+  a  ( l,  P )  ( "h  V  kprp  ^  ^d/d9 )  ( 1 p\( cose ) ))' s i  n<t> 

■i  pl(cose) 

■  i3(4p)('1/knrn)((3/3rHrh/k,r)))  1 


p  P' 


l  p  ' ' '  sine 
r=r_ 


]  , 


(56) 


<Ve 


V  0 


l  1 

£=1 


;l  21+1 

wm 


[b(£,p)(_ 


P^(cose) 
sine 


-)(sin^)j£(kprp) 


+  ia(^p)^1/kprp)((3/8r)(rj^kpr)))(d/d0)(P1t(cose))sin<ji 

pl(cose)  p 

+  PU,  p)( - iTRe-^3in^hJ(yp) 

+  ia(^p)(1/kprp)((3/3r)(rhi(kpr)))(d/de)(Pi(cosP))sin^ 

r=r 


(57) 


(Hp\  =  '  v  E°  J1i^^Wtb(^p)(“j£(kprp))((d/de)(pl(cos0)))cos* 

pl(cose) 

+  ,i'(£,p)(1/kprP)((!l/8r)<rJi(V)))r=r(J5Tm-  >cos* 

P 

+  3(>e,p)(_h£(kprp^  (  (d/d0 )  (P^(cose) )  )cos<t> 

,  pl(cose) 

+  i9(4p)(-1/Vp)((3/3r)^V))J=r(-W')cos^  <58> 
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The  boundary  conditions  implying  continuity  of  the  tangential 
component  of  the  electric  vector  in  the  0-direction  may  be  described 
by  the  rule 


IT 


oo 


(4p) 


P^(cose) 

sine 


iBU,p)((d/d9)Pi(cos6)) 


P»(cose)  , 

+  A(£,P)^W  "  i8Ujp)((d/de)P^(cos9))] 

“  pl(cose)  , 

=  ^0  Jj^p+D-sTne-’  iBU,Pn)((d/d9)P£(cose>> 

pl(cose)  1 

+  Vp+l)  sine  -  iS(£>p+1)((d/de)P£(cose))]. 


(59) 


Here 


r  (  f>)  -  <  2£+l 

CU)  1  ink)  * 

(60) 

flUP)  *  c(a^(kprp)a(4p>  • 

(61) 

nu,p+D  ■  c<*)jVVirp)au.p+o  • 

(62) 

A(-£,p)  ■  W)h\%r^u,p*l)  ' 

(63) 

AU  P+1)  =  CU)h£(Vlrp)'‘(«,p+l)  • 

(64) 
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BU,p)  =  C(i){l/kprp)(a/ar)trhJ(kpr>)  8(£>p)  ,  (65) 

r_rp 

Wi) =  c(£»1/Virp)(s/9r,(rhil(Vir)’Lr  6«.p«)  •  (66) 

P 

BU,P)  =  C^Hl/kprp)(3/3r)(rj£(kpr))^b(£jp)  ,  (67) 

P 

bu.p+d  =  c(£)(1/Virp)(3/3r){rj^(Vir^rb(^P+i)-  (68) 

p 

Letting 

SU,p)  =  AU,p+l)  +  A(£,p+1)  '  AU,p)  ‘  A(£,p)  (69) 

and 

Tu,p)  =  BU,p+i)  +  8U,p+i)  '  BU,p)  "  Bu, p)  ’  (70) 

we  deduce  that 

»  pl(cose)  ! 

jjISft.p)  -W-  -  1TU.p)(8/a,,lVcoso)1  -  0.  (71) 
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implies  that 


Observe  that  (EJ  =  (E  ) 

<i>P  vVp+1 


oo 

*Eo  ^1^"A(^,p))(d/d0)pi(cose)-i(-E 

+(-*kP))(d/de)p^cose)“1<-ButP)) 


+  ^P^(cose) 
U.p}^  sine 


P»(cose) 

JL^ - 1 

sine  J 


nEo  J1[(-AU,p-n))(d/de)pi(cose)'i(-B(£>p+i)) 


P^(cose) 

sine 


+(_^fp+l)^d/de>pi(cose)-l(-8^p+1)) 


P^(cose) 

sine 


-]  . 


(72) 


Thus, 


llt^£,p)'Wd6)^(C0Se)1'n(4p)l!^’> 


t 


and  we  conclude  that  S(^p)  =  o  and  T{£  p)  =  o. 


Upon  introducing  the  constants 


AUpf  a(4p)(1/kprp)(3/9r^rj^kpr))  C(^’ 

K  r=rp 

% 

b(4P)  =  b(4p)jvvyc<^* 


(73) 


(74) 

(75) 
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(76) 


^P)  =  a(^P)(1^pV)(3/3r)(rh]e(kPr^r  CU)’ 


2 


T 

U.p)  =  SU,p)h*(kprp)C(<>)  , 


Au,  p+i) =  a(£,p+i)^/{VirP))(3/3r)(r^(Vir)J=crU) 


(77) 

(78) 


b(4p+1)  =  bU,p+l)j£(kp+lV  C(£)' 


(79) 


l(^P+i) =  “(4P+i)(1Vp))(3/5r)(rh£(Vir)J=^) 


(80) 


Wl)  =  e(^P+l)h1£(kp+1rp)C(^), 


(81) 


and  setting  (Hp)0  =  (Hp+1)0,  we  arrive  at  the  equality 

-(^1*Eo  +  «Up) «d/-.)Pj(eo..)) 

,y.  Pp(cose)  -v,  i 

+  E»,p)<-  -W-}  +  w  (^pjttd/deJP^cose))] 

k  ,,  oo  -v  p^(cose)  v  i 

■  -  ^r"Eo  j,W-  +  iWi)Sd/de,p*(cose)) 


+  8(^P+1)( 


pl(cose)  n  i 

+  i'WnWd9)p<tc°s9))1- 


(82) 
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r 


Now  set 


sU,p)  =  ®(£,p+Dkp+i  +  ®U,P+I)kp+1  *  Bu.p)kp  '  Bu,p)kp  (83) 


and 


Tu,P) =  V,p+i)Vi +  ^u,P+i)kp+i "  A\t, P)kp  ■  Au,p)kP 


(84) 


Equations  82-84  yield 


«°  -v  pj(cose)  ~  , 

Jj[SU,p)(-  '  sine  )  +  'TU  p)(d/do)Pt(coSe)] 


=  0 


(85) 


The  boundary  condition 


‘W* 


(86) 


Is  now  utilized. 

We  observe  that 


,  k  «  <v.  ,  v  pj(cose) 

(l^)lrE0  It-B(£,p)((d/d9)P£(COS0))  +  iAU,p)(" 


sine 


,  pj(cose) 

'  +  MU.P)(-W-)I  *  <HpV 
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II  r — i  8 


00 


and  (Hp+1^  = 


“  pj(cose) 

*Eo  p+l)^d/de >p^(cos0>)  +  lA(£,p+i)^  sine  ) 


~  P:(cose) 

5(£,p+l)^d/de^P£^COs6^  +  1  (£,p+l)  sine  ^  ■ 


(87) 


Thus, 


00  f\j 

Jx[_  BU,P+i)kP+i 


U,p+I)kp+1  +  BU,P)kP 


8U,p)kp]‘d/d0)pl(cos0) 


+  i 


,n^,p+i)Vi +  Au,p+i)kP+i  -  V,p)kp 


■  V.p)kp]( 


^(cose) 
sine  ' 


=  0  . 


(88) 


This  implies  that 


_  1 


^  ^  pj(cose) 

[su>p)((d/de)P£(cose))  - 


(89) 


0/ 

Now  equations  85  and  89  are  used  to  determine  the  values  of  pj 

and  T^  Multiplying  both  sides  of  equation  85  by  P^i  (cose)slne 

and  both  sides  of  equation  89  by  [ ( d/de ) , (cose) ]sine  and  integrating 
from  0  to  it,  results  in  equation 
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00  ^ 

£=1 


f-n  pJ(cose)P»,  (cose)  ,  . 

[— - 5-^ -  +  ((d/de )P» (cose ) )  ( (d/de )pi ,  (cose ) )  [  si  node 

Jo  sin  e  'L 


-i^lT(£ip)  [  fP£(cose)((d/de)P^,(cose)) 

+  ((d/de)P^(cos9))P^,  (cose)3de  =  0 


(90) 


which  implies,  in  view  of  equations  51  and  52,  that  ^  =  0. 

symmetry  =  0.  We  conclude  that 


By 


B(f.Ptl)hf(lVlrp',lVl  *  b(£,p«-I)VVl'"plkP'U 

=  B«,p)h)(kprp)kp +  bu,p)¥kprp)kp  ■  <91> 


Now  the  associated  relation  derived  from  equating  the  tangential 
components  of  the  t  vector  is,  from  equations  71  and  73,  the 
following: 


1(s/3r>(rJ*Vir)”r>,pH) 


*  (o:)[(w,ri(rtl<kp.i'i;;r  ,p+n 
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■  iTF— i b(£,P) 
P  P  p 

+  (o-)'<8/sr)(rhi(kpr)>'  bu,p) 


Remark.  If 


then  the  inverse  of  A  is  given  hy 


A_1= 


a22/ A  -a^/ ^ 


\-a21/A  a^/Ai 


where  A  =  "  ^12^21  * 


We  define  for  the  sake  of  economy  several  terms,  following 
Shapiro  et  al .  (13).  Let 

Etf,p)  5  Cfc  ^/^MrhlCKpr))  =  rV-  ^<Vp)  *  WIVp» 


r=rp  P  P 


=  r\r-{8/3  P)(phj(p))_ . 

Vp  1  p  kprp  ’ 


{<‘.p+i>  *^k(a/ap)(ph‘(p)H^P 


u,p)  ■  (»/ap)(pj,(p))0,k 


P  P 


r 

P  P 


*  FT"  0/*-'>(rj|t(k  r)) 

P  P  H  P 


1 


(*,p+d=  {3/9p)(pj«(p))p= 


'p+i  p 


p'kP+irP 


J(*,p)  =  Wp>  • 


ju,p+l)  ^l^pvl'v  • 


hu,p)  "  "IVp’  • 


(96) 

(97) 

(98) 

(99) 

(100) 

(101) 


and 

h(5,  ,p)  ~  Vkp+Irp^ 


(102) 


Now  the  relations  between  the  coefficients  in  matrix  form  are 


(jU,p+l)kp  +  l 

nU,p-*-l) 


=  B(^P  +  1) 


CU>P+1)  J 


(103) 
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ju,p)kp 


'(£»p) 


hu,p)kp 


'U,  p) 


=  B(*,P)  ■ 


Observe  that 


u. p)  (£,p)|  =  ru,p+d  lb(l> p+i) 


Hi,  p) 


*U.P+1)- 


Thus,  letting 


RU,P)  ’  (b'£-p))-1b«-"+i)  , 


we  deduce  that  since 


(B(£,P))-1  = 


CU,P)/kpAp 


rnu.p)/kPAP 


■hU,P)/AP 


J’u,p)/Ap 


(104) 


(105) 


(106) 


(107) 


where 


(108) 


(109) 


r(£,p)  . 

r(i,d  ■ 

[!(i,p)JU,ptl)(-^i)  '  hW,p)n«,p+i)I/sp  • 

(110) 

dU.p)  _ 
(1,2)  " 

+  k  ,  + 

[htt,P+i)5tt,P)(k^  .  hU,p)^U,p+i)]/Ap  ’ 

(111) 

RU, p)  . 
(2,1)  ' 

+  -  k  , 

[-nU,P)jU,P+i)V  +  nU,p+i)j(£,p)^/Ap  , 

(112) 

d(*,P)  = 
'(2,2) 

k 

[-hu,p+i)nkp)("^i)  +  5u,P+i)jU,p)1/Ap  • 

(113) 

Now  on  to  the  development  of  the  matrices  relating  the  a-a 
coefficients  in  layer  p  to  those  in  layer  p+1.  The  first  relation 
is  derived  from  equation  69  and  the  definitions  60-68,  and  is  expressed, 
using  the  notation  in  equations  95-102,  as 


V,  p+1  fee,  p+1)  +  j(£,p+i)a(e,p) 


'  h(£,p)a(£,p) 


J(e,p)a(e,P)  =  o  . 


(114) 


The  next  relation,  derived  from  equations  74-81  and  equations  83-90, 
takes  the  form  of 


“(£, P+1)  kp+|rp(3/9r)(rhe(kp+lr)J=r  C(£)kp+1 


+  a(e,p+D  k^(3/9r)(rje(kp+ir)J=r_c(£)Vi 


■a^p)^(8/*r)(rhi(kpr)]LrcU)kp 

M  ^  p 


-  a 


(l, p)  r))^  CU)kp  =  0 


(115) 


which,  after  using  the  notation  expressed  by  equations  95-102,  may  be 
written  as 
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'U,p+l)^,p+l)kp+l  +  a(£,p+l)nU,p+l)kp+i 

■  V,p)cu,p)kp  ■  aU,p)nU,p)kp  =  0 


(116) 


Let  us  define 


iU,P) 


l3U.  P) 

hkp)  \ 

lV,P)kp 

cu,p)V 

and 


AU,P+1) 


Define 


h(£,P+l)  \ 

ln(£,p+l)kp+l  ?U,p+l)kp+l/ 


(117) 


(118) 


QU'P)  =  (A^'P))"1  (A^'P+1)) 


(119) 
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and 
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where 


=  t5(t,p)Jlt,P+l)  ■  (“^i)hu,p)’’u,p«)1/a(i,p)  • 

Q(i'J)  ■  "((.pfepti)  ■  (“^1),’kp)cK,pti)I/a«,p)  ’ 

H'.i]  =  [(“i^)ju,p)nu.p+i)  '  nU,p)j(-e,P+i)]/Au,p)  • 


(124) 

(125) 

(126) 


and 


n(£»p)  =  T(  C  *7 

Q(2,2)  kp  ,JU,p)5U.p+l) 


-  n 


U,p)h(£,p+1)]/A(£,p) 


(127) 


Now  we  wish  to  use  the  transition  matrices  and  R^,p^ 

to  get  relations  between  the  internal  and  the  external  coefficients. 
First,  note  that  a(£>1)  =  •  0  and  a^N)  =  b(£N)  =  1,  where 

N  is  the  number  of  regions  into  which  space  is  subdivided  by  N-l 
spheres.  Second,  observe  that 


lau,n 

1 

.  q(£,i)q(£»2)  ...  qU.N-D 

l  - 

L  0  J 

“U,  N) 

(128) 
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and 


’u,n 

o 


.rU.1)rU,2>  ...  RU."-1) 


U,n) 


(129) 


or  setting 


Q  =  q(^>  1)q(-^» 2)  ...  o(^'N-1) 


(130) 


and 


R  =  rW.1)rU.2)  ...  rU.N-D  , 


we  have  the  following  relations 

/« 

\  0  /  VQ 

and 


/ 

aU,l) 

0 


\.r\ 


(1,1) 

Q(l»2)\ 

!(2,1) 

Q(2,2 )/ 

'(1,1) 

R( 1 ,2 ) 

'(2,1) 

"(2,2)/ 

I  X 


\ 


(131) 


(132) 


(133) 
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Thus,  we  see  that 


att,n)  *  "Q(2,l)/Q(2 


.2) 


and 


(134) 


(41)  Q(l,l)  "  Q(1,2)Q(2,1)/Q(2,2)  * 


(135) 


Furthermore,  once  a 


“(4 P+1)  and  aU,p+i)  by  the  relation 


(4 p)  and  a(^,p)  are  determined,  we  obtain 


hu,p)\ 

_  ! 

\a(4P)y 

\ 

q(4p) 

g(l,l) 


i(4p) 


(136) 


Also,  we  deduce  from  equation  129  that 


SU,N)  *  R(2,i)/R(2i2) 

and 


(137) 


(4D  R(l»l)  ~  R(1,2)R(2,1)/R(2,2)  * 

As  before,  once  e^p)  and  ^  are  determined,  we  obtain 
8U,p+l)  and  b(^p+i)  by  the  relation 
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/bU,p)\ 

\V,p)/ 


p(^»  P) 1 

R(l»2) 


RU,P) 

”(2,2) 


\l,  P+l)\ 


fu,p+i)l 


(139) 


By  repeated  application  of  matrix  equations  136  and  139,  we  determine 
all  expansion  coefficients;  and  thus,  using  equations  53  and  54,  com- 
pletely  determine  electric  field  E  and  magnetic  field  H. 


Determination  of  Total  Absorbed  Power 

The  Poynting  vector  is  generally  interpreted  as  a  vector  having 
length  equal  to  the  power  per  unit  area  traveling  across  a  surface 
normal  to  the  vector  and  direction  of  the  power  flow.  One  can  show 
by  the  Gauss  divergence  theorem  that  the  Poynting  vector  is  given  by 

•*•*■*■*■*■ 

t  ExH+ExH.  (140) 

b  -  2 


Maxwell's  equations  then  imply  that 

dlv(S)  ♦  y?  *  jjsfegl)  =  0  ,  (141) 

Poynting* s  theorem  in  differential  form  in  the  absence  of  impressed 
electromotive  forces  for  linear  material  media  and  where  e,  a,  and 
J  are  the  permittivity,  permeability,  conductivity,  and  electric 
current  density,  respectively,  and  the  sign  *  attached  to  a  vector 
denotes  its  complex  conjugate. 

Let  us  write  (in  terms  of  the  spherical  coordinate  base  vectors 

V  V  and  V 
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E  =  E  e  +  EA  eA+  E  e 

r  r  6  6  <k  <t 


(142) 


and 


where 


H  -  Hrcr  +  H0ee  +  V*  - 


(143) 


er  -  si necos<{)i  +  s i nes i r»4> j  +  cosek  , 


e.  =  cosecos<(ii  +  cosesin<j>j  -  sinek  , 
u 


(144) 

(145) 


and 


e  =  -  sin<(>i  +  cos<}> j  . 


(146) 


Observe  that 


er  x  ee 


sinecosifi  sinesint  cose 

cosecosij!  cose  sin<t>  -sine 

2  2  2 
i (-sin  esin<j>  -  cos  esin<$>)  -  j(-sin  ecosij* 
-+ 

+  k(sinecosesintcosi}i  -  sinecosesin<Kos<j>)  = 


=  e,  . 


A  similar  calculation  shows  that 


ee  x  %  =  er 


-  cos^ecos<j>) 

-►  -> 

-si  n*}>i  +  cos<f>j 
(147) 

(148) 
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(150 


E  x  H  =  M  '  HeVer  +  <EA  ’  \*rK  +  <ErHe  “  HrEe>ee  • 


What  we  need  to  compute  is 


S  •  (-N)  =  (E  x  H)  •  (-N) 


(15i: 


When  N  =  er  .  Now  the  power  going  into  the  sphere  is 


r*  (2ir 

[ 

0 


-  (EeH*  -  He^)sinede]d^  . 


(152) 


At  this  point,  we  stop  to  refamiliarize  ourselves  with  the  struc¬ 
ture  of  the  vector  wave  functions  listed  below: 


:(ej) 

M(l,n) 


1 


Sine  ^V)Pi(cose)si,i*e8 


zj(kpr)(d/de)(pJ(cose))cos<f.e(f  , 


(153; 


M(l,n) 


^^(kpr)pJ(cose)cos^e 


z^(kpr)(d/d9)(pJ(cos9))sini»e^  , 


(154 


•+(o,j)  _  njn+l).  2J(k  r)pj;(cos9  )sin$e 
N(l,n)  kpr  nv  p  nK  r 

+  fV  (rz^(k  r))(d/de)(P^(c°se))sin<t>e6 

P 

+  ir-FiTF9<8/3r) (rzn(kpr) )Pp(cose )cos^4,  ’ 

P 

d 

N(l,’n)  =  zn(kpr)Pn(c0S9)c0S^r 

+  i^r  (3/3r)(rzj(kpr))(d/d9)(P*(cos9))  cos  <j>eQ 

-  Tn^3/,r><rz^y)|p^cosels1nn  ; 

p 


to  contemplate  the  use  of  the  fact  that  in  the  region  (p= 
the  biological  material  the  electromagnetic  field  is 


E  =  E  y  in  rM (° » 1 )  .  *w(e»l) 

ony  HTH+lJ  lM(l,n)  lN(l,n) 

H  =  -  \  r  V  in  2|3+1.  _  rM(e  , 1 )  +  iN (o >  1 ) 

w0«  Eo  ^  nTn+IJ  [M(l,n)  lN(l,n) 

+B(n.»)Sii:nj+,«(n,N)"j?:nil* 


(155) 


(156) 


surrounding 


(157) 


(158) 


or  more  compactly 


-+•  -+ 

E  =  E1  +  £r 


and 


H  =  H1  +  Hr 


(159) 


(160) 


and  in  accordance  with 
of  %  has  been  deleted 
overbar  — ) 


Stratton  (14,  p.  568),  we  write  (where  a  factor 
and  now  the  complex  conjugate  is  indicated  by 


S  = 
r 


E.H.  ' 

^  $  (J>  0 


(161) 


and  take 


W3  =  -Re 

a 


fir 


■2n 


-  2  . 

Sr  r  smed(j)]d9. 


(162) 


o  o 


Observe  that 

-*■  ±  +  -*•  -V  -V 

h(E  x  H  +  E  x  H)  =  Re(E  x  H) 

-*■ . 

r"1  ^  ni 


=  Re(E  x  H1)  +  Re(E1  x  Tf  +  Er  x  H1') 

— V 

+  Re(Er  x  \f). 


(163) 


A  direct  calculation  shows  that 

0  Re(E1  x  T?1’ )  .  (-N)dA  -  0  . 


(164) 
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Thus,  to  get  the  total  energy  absorbed  by  the  sphere,  we  need  only 
compute 

Wa  =  Wt  -  W$,  (165) 

where  represents  the  energy  dissipated  as  heat  plus  the  scattered 
energy,  and  Wg  represents  the  scattered  energy. 

Now 


Ws  =  Re  f*  fV*  "  EAHe)r2si‘n9d0d^ 

>  o  o  <*’  ^ 

2tr  7T  00 

" Re  L  LlE°  Jif£ 

’ Re L  L'e°  J/' i6u,N)<N!?:l)Hj 

•  '^Fo  jj5  sff^s.NlKM)^  -  ,“(s,N)(R(^)>41,r2si"8d* 


(2 ir  fir  °° 

Re  l 

•*0  O  £=1 

kw  «»  _ 

{-  — ^  r  1  is 

V  °  S=1 

( 2ir  fTT  00 


-  I  I 


(2S+1 ) (21+1 )  rc(g,g) 


^  ^  sts+i1tlMttF(^s)aU,N)a(s,N)  +  F(I,’sjaU,N)B(s,N) 


x  r(<*>£5)  t 

T  r  _  \<*/  »  H\P| 


+  p(3>Ct)o  “  +  c(B,BL  O  1 

FU,s)BU,N)a(s,N)  F(£,s)P(£,N)e(s,N)]  ’ 


(166) 
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Observe  that 


,2’  r'-!V 

U  0 


'  r>  ■> 


0  0 


(£,N)a(s,N) 


•  (j~ r  (3/3r)(rhJ(kNr»(d/de)P^(cose)sin(fr 


+  (j~  (3/3r)(rh^(kNr))(d/de)pJ(cose)cos<{l) 


•  (|y-sine  ^a/3 r)  (rhs  (kNr ) )ps (c°se  )cos$)Jsineded$ 


=  R  n  c(e»a) 

V,N)  (s,N)Ffc,s)  • 


Since  for  all  positive  integers  l  and  s,  we  have 


[P^(cose)(d/d9)P*(cose)  -  p*(cose)(d/d9)p|(cose)]de  =  0, 


it,  thus,  follows  from  equation  170  that 


c  (  ^  »®  )  -  n 

?U, S)  "  0  1 


(170) 


(171) 


(172) 
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and  an  almost  identical  argument  shows  that 


Ffes)  ■  0  <173) 

for  all  l  and  s.  Note  that  the  value  of  W$  is  independent  of 
r.  This  enables  us  to  make  use  of  the  asymptotic  formulas 


hj(p)  4  (-try* 


(174) 


and 


(d/dp)hj(p)  «  -  -y(-i }n+1e1p  +  ^-i)n+1eip  .  (175) 

P 


Further,  we  have 


p(a»a)  _  f2  f  ^2  -S+£/  .  ,S  r/u(o,3) »  (_■!  Uw(°»3) ) 

FU.S)  J0  Jo  v  Eo  ’  <  »  (  1,(N(l.s)>( 


'(M(i!J))6('i)(N(°’='VrZs1neded* 


ri t  irk 


J-?E l  iS+£(-l)S(-i)[(M^;|j)  (N[^)) 


u  to  0 
o  o 


ir2s("sde ; 
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r(a,a) 

FU,s) 


=  \]  iS+£(-l)S(-i)[h^(kNr)(l/kNr)(3/3r)(rh|(kNr)) 

1  1  Pp (cOS6 )P^ (COS0 )  , 

-( (d/de )P^(cose ) ) ( (d/de )P*(cose ) )  -  — - —I - ]r2 


.  ~T 

sin  e 


sinede 


!^j  r2  .-s+£+l/  ,xS  2  (s+1 ) !  ,  ,1U 

vou  Eo  1  ^  ^  2s+l  (s-1)!  S^S  ^6(s,£) 


*  ^(kNr)  (l/kNr)  (3/3  r)  (rh*(kNr) ) 


(176) 


after  applying  Proposition  6  for  m=l. 

To  complete  the  calculation,  we  make  use  of  the  following  lemma. 

Lemma  1 .  For  all  k , 

h^(kr)(l/kr)(3/3r)(rh*(kr))r2  =  -i/k2.  (177) 

Proof.  Equation  177  is  equivalent  to 


MkrH£  h£(k>")  +  r^h£'(kr)]  =  -i/k 2  .  (178) 

In  view  of  equations  174  and  175,  this  completes  the  proof  of  Lemma  1. 
Using  Lemma  1  and  equation  176  we  deduce  that 

i i2s(-l)si (-i )2s2(s+l)2 


\V,l\ 

_  /kN, 

=  (— ^] 
V 

2 

E  5 to  - \ l  V ■ 

-0-(.4sJ.  1 

K2(2s+1) 
■2.  .,,2, 


(179) 
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Similarly 


F(il!  ■  <lnr>2Eo  V,s)s2(s+1)2/'(2s',1)kN)  . 


(180) 


=  OE0  j  <2stl><l»(s,N)l2  +  lS(s,N)'2>/kN 


(181) 


Energy  balance  is  maintained  through  the  introduction  of  a  third 
term,  namely 

"t  ■  -Re  f*  f  (Ee  *  Ee  "I  -  E[  "e  -  EX>s<"8«M«  .  (182) 

0  >  0 

Thus,  collecting  multiples  of  the  expansion  coefficients,  we  obtain 

wt  ■  -Re(f  1’ '  j/  zfray  [“(i,N)<Mi?;l!>-  1  WImI’,1 


•  jj'-1’5  iH^rTl(M»:s)V  i(No:s)V 
+  j/ zf&iy«Ma:l)>e  -  1("f?:l)V 
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where 


FU,s) 


-  1)dA 


(0,1). 


(184) 
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t 


c(l.o) 

FU,s) 


-  km[;$)  -  KNj*’}i)  H-i)(N }dA, 


'(I,*) 


(1,S)' 


(185) 


F«.’sr  L“«iw}>  -  *<!!>  i 


(186) 


-  <"(*:!)>  w(WJ)  -  ’wI&Jmm*- 


(187) 


First  we  compute  Observe  that  if  we  let 


Afes)  ■  Io  10  IWh|(kNr)Pl(cos9)cos* 

[-  jc(kNr){(d/de)pJ(cose))cos4> 


1  1  o 

-  kj^VsTne (3/3r)  (ns(kN'r) )P$(cose )cos<|>]r  sineded* 


AU;s)  -  ^(kNr)(-js(kNr)) 


rlT 


P^cose ) ( (d/de )P^ (cose ) )de 


„  ,  -  t*  p}(cose)P*(cose) 

-  i1tr2h^(kNr)(1^F)(3/3r)(rjs(kNr))  J  - ^ - 


and 


B 


(cl)  _ 

U,  S) 


2ir  n 


-  hj( kNr)(( d/de  )pJ (cose  ))sin(j>[ 


-  ?lW  Js^Nr^ps^cos0^sin* 


-  (3/3r) (rj$(kNr)  )((d/do )P^ (coso ))sin<}>] r^sinedo 


2  . 


r 


((d/de  )pJ  (cose  ))P*  (cose  )do 


then 


■  "■2h1t(kNr)js(V) 

+  hJ(kNr)(a/3r)(rjs(kNr)) 
.  (d/de)P^(cose)]sinede, 


r(^t»l)  _  n  ( (*  j  1  )  d  ( '-*  >  1  ) 

F(P.,s)  -  Au,s)  '  B(f  ,s) 


In  view  of  equations  18  and  170,  we  obtain 


[(d/de)pJ(cose) 


1 


de  (188) 


(189) 


(190) 
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F 


(a.l) 

U,  S) 


(a.l) 

U.  S) 


B 


(a,l) 

U,s) 


=  -*r2h£(kNr)js(kNr) 


r  it 


[P,(cose)((d/de)P*(cose)) 


+  P^(cose)((d/d6)P^(cose)3de 


iirr2h^(kNr)(3/3r)(rjs(kNr))  p  P^(cose)P^(cose) 

J  ft 


.  2 
sin  e 


+  ( (d/d0 )P^(cose ) ) ( (d/de )P* (cose ) )]sinede 


1 7cr^h  ^(  kNr  )(3/3r)(rjs(k|«jr)) 


kNr 


\1,S) 


2s+l  Ts-lT! 


s+1) ! 


s(s+l) 


iirrh^(kNr)(3/3r)(rjs(kNr)) 


6  2s2(s+l)2 

6 {l, s)  2s+l 


(191) 


Next  we  compute  F 


(l,a) 
U,  s) 


.  As  before,  we  let 


(192) 


and 


B(l»a)  _  .• 

B«.S)  '  •’ 


(193) 
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Use  of  equations  153-156  yields 


V,s)  ' 


2t r 


j  t?WJ£(kH'-|Pl(C0Se>C0S't 

'  0 


(3/8t”)(rj^(kNr))(d/de)pJ(cos0)cos<^] 


•  (a/3r)(rhJ(kNr))Pg(cos6)cos<jir2sineded<j> 


=  -  ^  j/(kNr)(a/3r)(rhJ(kNr))  [  p£(cose)Ps(cose) 
N  -’o  sine 


de 


1  7T 


Y  ((3/3r)(r3£(kNr)))((3/3r)(rh^(kNr))) 


( (d/de )P^(cose ) )P* (cose )de 


Also,  using  equations  153-156  and  193,  we  deduce  that 


>(l»a) 

*(*.  S) 


/•2tt  rit 


D/ n  ~  \  I 


r  i 

[  ( ~  J^(  kN»") )  (d/de  )P^(cose )  si  n<t> 


+  kjjjrsine  (3/30(n£(kNr))P*(cose)sin<j>) 


•  (3/3r)(rhg(kNr))(d/de)P^(cose)sin4>]r2sineded<f> 


(194) 
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Carrying  out  the  integration  with  respect  to  <f>  we  see  that 


Bfe)  ■  if  J£0‘Nr)C3/3r)(rh|<kNr)) 

•  |  ( (d/de )P^(cose) )((d/de )P^ (cose ) )si nede 

+  ^  ((3/3r)(rj£(kNr)) )(3/3r)(rh^(kNr))) 
kN 

r  i  i 

•  Pi (cose ) (d/de )P^(cose )de . 

Jo  *-  s 

From  equations  185,  194,  and  195,  it  follows  that 


(195) 


F(l,a) 

Vs) 


A(i,°0 

Vs) 


R(l»a) 

Vs) 


Thus , 


(196) 


Ffe)  =  '  VkNr)(8/ar>(rh!(kNr>> 

fir  p}(cose )P*(cose )  ,  i 

•  [— - ^ - +  ((d/de  )Pi(cose))((d/de)P^(cose))]  si  nede 

o  sin^e  c 


-  ^((3/3r)(rj£(kNr)))((3/3r)(rh^(kNr))) 
kN 

•  t((d/de)(p)(cose)))P1(cose)  +  ((d/de)P*(cose))p](cose)]de 
J0  c  s 

=  "  kNr ^ 3 /3 r T'hs ( kNr } ) 16 (^, s ) 2|s+!'+1  ]  •  (197> 


Furthermore 


Re(a,  ^&A\  +  a,.  k,sF^*“b 


l(s,N)  (s,s)  T  (s,N)  (s,s) 


=  Re  [a 


(s,N) 


-irrhg(kNr)(8/3r)(rjs(kNr)) 


2s2(s+l)2 


'U,  s)  2s+l 


+  a 


(s,N) 


-inrjs(kNr)(3/3r)(rh^(kwr)) 


s^Nrjj  .  2s2(s+1)2  , 

-  6(l, s)  2s^l  "  J 


,  ,  2s2(s+1)2 

"  Re<a(s,N)^  2s+l  ' 


Hence,  using  the  fact  that  k^  is  real,  we  have 
2ttE2  , _ 

“t  ■  -T-VTaT  Be  I  (2s+l)(o(s^)  ♦  S(s  N)) 

K»,  S  I 


This  follows  from  an  induction  argument  and  the  fact  that  if 


/  ,-\n+l 

u  +  iv  =  A — f- —  [cos(kNrl  +  i  sin(kNr)] 


*Vi 


and 


w  =  y~  cos[k^r  -  , 


then  for  all  real  numbers  A  and  B, 


(198) 


(199) 


(200) 


(201) 


1 


1 


1 
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Re{  i  [w1  (u+iv)(A+iB)  +  w(u  '-iv*)(A-iB)]> 


=  (v'w-vw')A+  (u'w-uw')B 


A_ 

k. 


(202) 


for  every  positive  integer  n.  A  prime  on  u,  v,  or  w  denotes 
differentiation  with  respect  to  r. 

Thus,  time  averaging  shows  that  the  total  absorbed  power  is  given  by 


W 


a 


It°V^0  V2n+1)(a(n,N)  +  8(n,N))! 
kN 


nE, 


II  '  *  -k 


I  ). 


(203) 


Summary  of  Key  Equations  and  Formulas 


In  summarizing,  we  set  down  the  key  equations  and  formulas  upon 
which  program  CSM  is  based. 

Fields  for  the  p-th  region: 


E 


P 


E  exp(-^t)  l 
0  1=  1 


:z  2.e+i 

1  KZ+T) 


(°a) 

U,p)m(i,-0 


[a to  _ \M 


ib 


u,p)r 


(e,l) 

(1,-0 


au,P)M(i’j)  •  i8u,P)N(i;')L 


(204) 
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B»  •  -  A  E„exp(-»t)  IgJjy 


P  VQW  ”0 


+  t,(e,3)  .  tj(0,3). 


(205) 


where  the  vector  wave  functions  j *  and 

are  obtained  by  replacing  the  spherical  Bessel  function  Jn ( ^pr )  by 
the  spherical  Hankel  function  hn^(kpr)  in  the  expressions  for  the 

vector  wave  functions  and  N(i'l)  ' 

Complex  propagation  constant  for  the  p-th  region: 

kp  *  Re(kp)  +  ilm(kp),  (206) 

where 

Re(k  J  =  |  [(1+  -i— ; 2fe)  )\  1]} 2  ,  (207) 

(%“)  P 

Im(kp)  =  |  t-f  l(l+  — ^ ,  (208) 

(e0<*>r  P 

-12 

Eq  *  free-space  permittivity;  =  8.85  x  10  F/m, 

ep  =  relative  dielectric  constant  of  p-th  region;  *  1 
for  free  space, 

Op  *  conductivity  of  the  p-th  region;  =  0  for  free  space, 
w  =  angular  frequency;  =  2it  x  frequency  (in  MHz), 

8 

c  =  velocity  of  light  in  free  space;  =  2.9979  x  10  m  . 
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The  field  expansion  coefficients  for  region  one,  inner  core  sphere, 
and  those  for  the  surrounding  medium  are  obtained  through  the  solution 
of  two  systems  of  equations.  Utilizing  the  notation  of  Shapiro  et  al. 
(13),  we  have,  with  ?  '<1  “(1,1)  "  eM)  "  »• 


aCt,  lj 

0 

11 

r-  — 

oij 

g-e,T 

—  - 

1 

1) 

_  o  _ 

It 

~  1 
R£,T 

'  1  " 

(209) 


(210) 


where  the  product  matrices  [ jl  and  [R^]  have  the  representation 


jiJ 

■f,T 


N-l 

n 

P»1 


0(4p) 

Q(iJ) 


(211) 


t,  T 


N-l 

n 

P=1 


d(-£>P) 

R(i  J) 


(212) 


with  each  factor  matrix  and  having  its  (i,j) 

elements  computed  by  means  of  the  following  formulas: 

Q(M) =  (A(f,p)rl[cu,p)ju,P+i)  '  ~k7hUp)n(f,p+i)]  •  (213) 


Q(1  >2)  =  {A(f,p))’lu(f,p)hU,p+l)  '  ^T^Up^Up+l)1  ’  (214) 
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(215) 


n(-t,p)  -  i  a  \-lf  P+1  j+  +  j-  i 

Q(2,l)  "  A(  £,p)  [  kp  J(^p)n(4p+1)  "  n(£,p)J(^p+l)]  * 

Q(2^2)  =  {A(£,p)}  1[_k^  j(4p)^(4p+l)  "  nUp)h(4p+l)]  ’  (216) 

R(M)  =  (A(£,p))1[_k^  c(4p)jU.p+i)  '  hU,p)nU,p+D1  •  (217) 

r(i>2)  =  (A(4p))"1[^IcU,p)hUp+i)  -  hUp)c(4p+l)]  *  (218) 

R(2j|  =  (A(^p)rl[jUp)n(4p+l)  '  “kj^  nU,p)J"(4P+l)]  ’  (219) 

R(2%P)  =  (A(£,p))’llJ'Up)C(4p+l)  '  “kf  nUp)hU,P+l)]  •  (220) 

Here 

A(^P)  =  JU,p)CU,p)  ‘  h(4p)nU,p)  ’  (221) 

jU,p)  =  J^kpr)  ’  (222) 

»Up)  =  Mkpr)  ’  (223) 

n(£,p)  =  2FT  t^1^(£-l,p)  "  %+I.f)1  ’  <224) 

CU,P)  =  WT  '  ^+1,^  ’  ^225^ 
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the  superscript  1  has  been  dropped  from  the  spherical  Hankel  functions 
h^(kpr)  =  jn(kpr)  +  iyn(kpr),  and  r  is  the  radius  of  the  boundary 
surface  of  the  p-th  region  for  subscripts  p  and  p+1. 

The  matrix  equations 


‘It,  if 

_ 

m%p+tf 

,  (226) 

1 — 

-%p+l> 

- 1 

fF 

n  (  P  ) 
K  /  •  •  \ 

(227) 

\  i  9  J  ) 

+ 

CL 

c or 

_ » 

» 

yield  the  expansion  coefficients  for  the  regions  p  =  2,...,N-1  in 
a  recursive  manner,  starting  with  derived  values  of  a^  ^  and  b^  ^ 
and  known  values  of  ^  and  ^  as  elements  in  the  left-hand 
members  of  the  matrix  equations,  and  employing  equations  213-225  for 
computing  the  necessary  coefficient  matrics  and 

Absorbed- power  density  at  an  interior  point  of  the  p-th  region: 


p  ■  o-VVV 


(228) 


where 

-»■ 

Ep  =  electric  vector  at  an  interior  point  of  the  p-th 
region. 

Op  =  conductivity  of  the  p-th  region, 

*  =  complex  conjugate  indicator. 

Average  absorbed-power  density: 

Pavg  '  <3'8’>'to/l'o)!S(EoVrN-l>'  <J29> 
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where 


Qa  ‘  lfi-Re  j  (WKxk*  'K.tdll-  !l  I(2m)(|a(t.N)lZ 

ku  -t-1  n  *- 

+  |8£>n!2)*  Qt  -  Qs  .  (230) 

eQ,Mo  -  free-space  permittivity  and  permeability  » 
k^  =  propagation  constant  of  the  surrounding  medium, 

=  scatter"'n9  coefficients. 


Total  absorbed  power: 

2P.  «>  2  2 

Pt  t  » -j1  [  (2t+l)[|lte(^  $  *  ^)l"  doj.^1  +  )1  • 

“  W  ’  (231) 

where  p2 

2 

=  power  incident  upon  a;  =(^-)rrr^_ i 


n  =  intrinsic  impedance  for  free  space;  =  376.7  ohms, 

rN  j  =  radius  of  spherical  surface  adjacent  to  the  surrounding 
medium, 

a  =  geometrical  cross  section  of  the  sphere  of  radius 
rN-l  *  =  27,rN-l/x’ 

A  =  wavelength  of  the  incident  wave. 


To  complete  our  summarization,  consideration  of  the  formulas  used 
in  generating  the  values  of  certain  functions  seems  appropriate.  The 
formulas 


pJ+1(cos8)  =  ^  cose  pj(cose)  -  ^  pj.^cose) 


(232) 


sine  (d/de  )Pj!j  (cose)  =  ncosepj(cose)  -  (n+DPj.^cose)  , 


(233) 
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together  with 


pj(cose)  =  sine. 

(234) 

P^cose)  =  3  sinecose 

(235) 

are  used  to  generate  function  and  derivative  values  of 
Legendre  functions. 

the  associated 

Special  limit  values  are  also  obtained  by 

Lim  Pn<cos9>  _  n(n+l) 
e-K>  sine  2  ’ 

(236) 

Lin,  Pn<cos9>  _  (-Dn+1n(n+l> 
e-*ir  sine  2 

(237) 

The  forward  recurrence  relation 

W2> +  vi<z>  ■  2n:‘  *„<2> 

(238) 

is  used  together  with  relations 

#  ,  cosz 

y0(z)  a  -  z  • 

(239) 

,_N  cosz  sinz 

n'2’  *  •  z2  -  z 

(240) 

to  generate  values  of  the  spherical  Neuman  functions. 

The  generating 

process  is  terminated  at  order  N  when  the  following  termination 
criterion 


ynU)|  i  STOPR 


(24.') 
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is  met.  Here  STOPR  is  a  number,  say  1.0D15.  The  user's  needs  will 
determine  whether  or  not  STOPR  should  retain  its  presently  suggested 
value.  Our  own  demands  were  satisfactorily  met  for  complex  argument, 
epr»  the  spherical  Neumann  functions  for  parameter  ranges: 

1.5  <_  |e  |  <_  1390.0  and  0.1  <_  r  <_  10  cm. 

The  backward  recurrence  relation 

Vi(z)  =ir1  Vz)  -  Wz)  (242) 

in  combination  with  an  appropriate  starting  value  is  used  to  generate 

values  of  the  spherical  Bessel  functions  of  the  first  kind,  i  (z) 

n 

This  technique  of  using  the  backward  relation  in  place  of  the  forward 
relation  helps  to  avoid  stability  problems. 

PROGRAM  DESCRIPTION 

Written  in  standard  FORTRAN  IV  for  the  IBM  360/65  system,  the  Con¬ 
centric  Spherical  Model  (CSM)  is  designed  to  calculate  the  internal 
absorbed-power  density  distributions,  average  absorbed-power  density, 
and  total  absorbed  power  for  a  spherical  shell  configuration  (simulat¬ 
ing  the  human  head)  subjected  to  plane-wave,  nonionizing  electromagnetic 
radiation.  Five  spherical  shells  plus  a  brain  core  sphere  are  generally 
treated,  but  provision  is  made  to  allow  as  many  as  eight  concentric 
shells  to  be  analyzed.  The  structural  components  of  the  head  model  are 
identified  by  regional  designators.  Regions  1  through  6  represent  the 
brain,  cerebrospinal  fluid  (CSF),  dura,  bone,  subcutaneous  fat,  and 
skin,  respectively.  Current  plotting  needs  of  the  USAF  School  of  Aero¬ 
space  Medicine  are  met  through  the  use  of  BGNSTP--an  advanced  general 
plotting  subroutine  package--and  the  use  of  the  CalComp  Model  936  digital 
incremental  plotter.  Since  the  plotting  package  is  not  available  for 
distribution,  the  plotter  calls  and  working  arrays  are  not  included  in 
this  published  version  of  the  program. 
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Basically,  program  CSM  consists  of  a  driver  routine,  five  sub¬ 
routine  subprograms,  and  one  function  subprogram.  These  routines  are 
single-entry  programs,  free  of  any  special  machine  dependence,  and 
utilize  subprograms  found  in  any  elementary  software  library.  Types 
REAL*8  and  C0MPLEX*16  signify  double-precision  real  and  complex 
variables,  respectively,  and  literal  data  appearing  in  a  FORMAT  state¬ 
ment  are  enclosed  in  apostrophes.  The  arithmetical  processes  are  per¬ 
formed  in  double-precision,  floating-point  mode.  This  feature  provides 
approximately  16.8  decimal  digits  of  precision  and  numbers  with  an 
exponent  range  of  -78  to  +75.  A  list  of  the  driver  routine  and  sub¬ 
programs,  including  function,  calling  sequence,  and  calling  arguments 
of  each  member,  follows. 

Driver  routine: 

Routine  MAIN  is  used  to  input/output  data;  to  compute  complex 
propagation  constants;  to  complete  the  calculations  for  the  absorbed- 
power  density  distributions,  average  absorbed-power  density,  and  total 
absorbed  power;  to  control  the  printing  activities;  and  to  direct  the 
course  of  calculations. 

Subroutine  subprograms; 

Subroutine  COEF  generates  the  expansion  coefficients  for  the  com- 
ponents  of  the  electric-field  vectors  Ep,  p  =  1,...,N0REG. 

The  calling  sequence  of  this  subroutine  is 
C0EF(ANP,  BNP,  ALPNP,  BETNP,  NMIN) 
where  the  calling  arguments  are 


ANP 


=  array  of  coefficients  for  vector  functions 


;(o,d 

( l»n) ! 


BNP 


=  array  of  coefficients  for  vector  functions 


(e.l) 

(1  ,n)  ’ 
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ALPNP  =  array  of  coefficients  for  vector  functions 
:(o,3) 


NMIN 


M 


(1 ,n)’ 


BETNP  =  array  of  coefficients  for  vector  functions 


N 


(e,3) 
(1 »n)  ’ 


=  number  of  terms  in  the  series  expansion  of 
each  component  of  the  electric- field  vector, 
E_  . 


The  above  arrays  are  double-precision,  complex,  and  each  array  is 
dimensioned  at  1000. 

Subroutine  EVEC  computes  the  radial, colatitude,  and  azimuthal 
components  and  the  scalar  product  Ep*Ep  the  electric-field  vectors 
Ep,p  =  1,...,N0REG. 

The  calling  sequence  of  this  subroutine  is 
EVEC(NP.PD) 
where  the  calling  arguments  are 

NP  =  region  identifier, 

PD  =  double-precision,  comp! ex, semi  completed 

absorbed-power  density  at  an  internal  point 
of  the  p-th  region. 

Subroutine  TERM  computes  (-l)n  or  (-l)n+1  times  the  appropriate 
part  of  the  n-th  term  in  the  series  expansion  of  each  component  of  the 
electric-field  vectors  Ep,p  =  1,...,N0REG. 

The  calling  sequence  of  this  subroutine  is 

TERM(NCK,  T,  KEY) 
where  the  calling  arguments  are 
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NCK  =  a  counter  count. 


T  =  part  of  the  n-th  term  in  the  series  expansion 
that  is  multiplied  by  the  appropriate  power 
of  -1  > 

KEY  =  0  for  T  to  be  multiplied  by  (-l)n  and  1 
for  T  to  be  multiplied  by  (-l)n+*. 

The  array  T  is  double-precision,  complex. 

Subroutine  BJYH  generates  the  spherical  Bessel  functions  jn(kpr), 
spherical  Neumann  functions  y^k^r),  and  spherical  Hankel  functions 

h<"(kpr). 

The  calling  sequence  of  this  subroutine  is 

BJYH(BJNP,  BHNP,  Z,  NN,  STOPR) 
where  the  calling  arguments  are 


BJNP  =  array  of  spherical  Bessel  functions  for 
p-th  region, 

BHNP  =  array  of  spherical  Hankel  functions  for  the 
p-th  region, 

Z  =  product  of  complex  propagation  constant  and 
radius  of  an  internal  point  or  boundary  sur¬ 
face  of  the  p-th  region, 

NN  =  maximum  order  of  the  spherical  functions, 

STOPR  =  a  test  quantity  for  terminating  the  gener¬ 
ation  of  the  spherical  Neumann  functions. 

The  arrays  BJNP  and  BHNP  are  double-precision,  complex,  and 
each  array  is  dimensioned  at  100.  Variable  Z  is  double-precision, 
complex. 


76 


Subroutine  PL  generates  the  associated  Legendre  functions  p^(cose) 
and  their  first  derivatives  with  respect  to  0. 

The  calling  sequence  for  this  subroutine  is 

PL (THETA,  N,  P,  DP) 
where  the  calling  arguments  are 

THETA  =  value  of  the  colatitude  angle  expressed 
in  radians, 

N  =  number  of  associated  Legendre  functions 

to  be  generated,  starting  with  the  func¬ 
tion  of  degree  one, 

P  =  array  of  values  of  the  associated  Legendre 

functions , 

DP  =  array  of  values  of  the  first  derivative 
of  the  associated  Legendre  functions. 

The  arrays  P  and  DP  are  double-precision,  real  and  are  dimen¬ 
sioned  at  101  and  100,  respectively.  THETA  is  a  double-precision, 
real  variable. 

Function  subprogram  MINN  determines  the  minimum  value  of  a  given 
array  of  positive  integers. 

The  calling  sequence  for  this  function  subprogram  is 
MINN(NRAY,N) 

where  the  calling  arguments  are 

NRAY  =  array  of  positive  integers, 

N  =  number  of  integers. 


The  array  NRAY  is  single-precision,  integer,  and  dimensioned  at 

10. 

Blank  COMMON  is  used  by  the  driver  routine,  MAIN,  and  the  sub¬ 
routines  COEF  and  EVEC.  The  list  of  the  arrays  and  variables  stored  in 
this  area  is 

★ 

FKP  =  wave  propagation  constants,  kp  , 

BJNP  =  spherical  Bessel  functions,  Jn(kpr)» 

BHNP*  =  spherical  Hankel  functions,  h|^(k  r), 

CEX*  =  exponential  value,  exp(-iwt),  for  circular 
frequency  u>  and  time  t, 

BDP  =  spherical  surface  boundaries, 

P  =  associated  Legendre  functions,  P^(cose), 

DP  =  first  derivative  of  the  associated  Legendre 

functions,  P^(cose), 

SIGP  =  conductivities,  o  , 

E0  =  intensity  of  the  incident  electric  field,  E, 

TIME  =  time, 

R  =  radius  of  an  internal  point  of  the  brain  core 

sphere,  or  spherical  shell  region,  or  the  radius 
of  a  spherical  boundary  surface, 

THETA  =  colatitude  angle,  e, 

PHI  =  azimuthal  angle,  4>, 

STOPR  =  a  test  quantity  for  terminating  the  generation 
of  the  spherical  Neumann  functions,  yn(kpr), 

NC  =  maximum  order  of  the  spherical  functions  minus 

2, 
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The  array  NRAY  is  single-precision,  integer,  and  dimensioned  at 

10. 

Blank  COMMON  is  used  by  the  driver  routine,  MAIN,  and  the  sub¬ 
routines  COEF  and  EVEC.  The  list  of  the  arrays  and  variables  stored  in 
this  area  is 

★ 

FKP  =  wave  propagation  constants,  kp  , 

BJNP  =  spherical  Bessel  functions,  Jn(kpr), 

BHNP*  =  spherical  Hankel  functions,  h^(kpr), 

CEX*  =  exponential  value,  exp(-iwt),  for  circular 
frequency  to  and  time  t, 

BDP  =  spherical  surface  boundaries, 

P  =  associated  Legendre  functions,  P^(cose), 

DP  =  first  derivative  of  the  associated  Legendre 

functions,  ^-P^(cose), 

SIGP  =  conductivities,  o  , 

E0  =  intensity  of  the  incident  electric  field,  E, 

TIME  =  time, 

R  =  radius  of  an  internal  point  of  the  brain  core 

sphere,  or  spherical  shell  region,  or  the  radius 
of  a  spherical  boundary  surface, 

THETA  =  colatitude  angle,  e, 

PHI  =  azimuthal  angle,  <p, 

STOPR  =  a  test  quantity  for  terminating  the  generation 
of  the  spherical  Neumann  functions,  yn(kpr), 

NC  =  maximum  order  of  the  spherical  functions  minus 

2, 
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NOREG  =  number  of  regions  in  the  Concentric  Spherical 
Model, 

NMIN  =  number  of  terms  in  the  series  expansions  of  the 

•V 

components  of  the  electric-field  vector,  . 

The  double-precision,  complex  arrays  and  variables  are  flagged  with 
an  asterisk  (*);  while  the  unflagged  arrays  and  variables  are  double¬ 
precision,  real— with  the  exception  of  the  last  three  members,  of  type 
INTEGER,  which  are  single-precision  variables. 

A  single-labeled  common  area,  COEF,  is  used  by  the  driver  routine, 
MAIN,  and  the  subroutine  EVEC,  for  values  of  the  expansion  coefficients 

a(n,p) ’  b;n,p)’  w(n,p) *  e(n,p)  stored  in  the  arrays  ANP’  BNP’  ALPNP’ 
and  BETNP  respectively. 

In  subroutine  BJYH,  if  variable  M,  the  maximum  order  of  the  spher¬ 
ical  Neumann  functions  y^(kpr)  (complex  kpr),  tests  <  2,  the  error 

message 


PROCESS  CAN  NOT  PROCEED  SINCE  NM2  =  0  FOR  Z  =  ... 
is  printed  out  and  the  computer  run  is  terminated. 

Input  to  program  CSM  is  by  keypunched  cards.  There  are  two  basic 
input  cards  with  structure  and  sequential  order  as  follows: 


Card  No.  1  (control  parameters) 

Columns:  1-10  FREQ.  Frequency  in  MHz.  (E10.3) 

11-20  EO.  Intensity  (field  strength)  of  the  in¬ 
cident  electric  field  in  volt/meter.  (E10.3) 
21-30  TIME.  Time  in  seconds.  (E10.3) 

31-40  STOPR.  A  test  quantity  for  terminating  the 

generation  of  the  spherical  Neumann  functions. 
A  suggested  value  is  1.0E15.  (E10.3) 

41-45  NORG.  Number  of  regions  in  the  concentric 
spherical  model  of  the  human  or  animal  head. 
(15) 

46-50  NOCR.  Number  of  cases.  (15) 
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Card  No.  2  (electrical  property  data) 

Columns:  1-10  EPSP(l).  Relative  dielectric  constant  for 

region  1.  (E10.3) 

11-20  SIGP(l).  Conductivity  for  region  1  in 
mho/meter.  (E10.3) 

21-30  EPSP (2 ) .  Relative  dielectric  constant  for 
region  2.  (E10.3) 

31-40  SIGP(2).  Conductivity  for  region  2  in 
mho/meter.  (E10.,3) 

41-50  . 

51-60  . 

61-70  EPSP(4).  Relative  dielectric  constant  for 
region  4.  (E10.3) 

71-80  SIGP(4) .  Conductivity  for  region  4  in 
mho/meter.  (E10.3) 

Card  3  is  a  similarly  structured  card  for  the  electrical  properties  of 
regions  5  and  6. 

Card  No.  4  (surface  boundary  data) 

Columns:  1-10  SBDP(l).  Radius  of  the  spherical  surface  for 

region  1  in  centimeters.  (E10.3) 

11-20  SBDP(2).  Radius  of  the  spherical  surface  for 
region  2  in  centimeters.  (E10.3) 

21-30  . 

31-40  . 

41-50  . 

51-60  SBDP(6).  Radius  of  the  spherical  surface  for 
region  6  in  centimeters.  (E10.3) 

Card  Nos.  5-(N0CR+4)  (coordinate  data) 

Columns:  1-5  NREG.  Region  number.  (15) 

6-15  R.  Radial  snnerical  coordinate  of  an  interior 
point  of  reqion  NREG  in  centimeters. 

Range:  0  <  R  <  SBDP(6).  (E10.3) 
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16-25  THETAD.  Colatitude  spherical  coordinate  of 
an  interior  point  of  region  NREG  in  degrees. 
Range:  0  <  THETAD  <  180.  (E10.3) 

26-35  PHID.  Azimuthal  spherical  coordinate  of  an 
interior  point  of  region  NREG  in  degrees. 
Range:  0  <  PHID  <  360.  (E10.3) 

The  last  card  of  a  single  data  set  must  be  a  termination  card  with 
the  symbols  /*  punched  in  columns  1  and  2.  Also  program  CSM  can  handle 
multiple  data  sets.  Each  data  set  [Cards  1 - ( N0CR+4 ) ]  is  stacked  one 
behind  the  other,  with  the  last  card  in  the  complete  data  deck  a  ter¬ 
mination  card. 

Program  printouts  consist  of  the  title 

ELECTROMAGNETIC  ENERGY  DEP0STI0N  IN  A  CONCENTRIC  SPHERICAL  MODEL 
OF  THE  HUMAN  OR  ANIMAL  HEAD 
followed  by  such  information  as 

FREQUENCY  =  ....  MHZ 
FIELD  STRENGTH  =  .  .  .  V/M 
TIME  =  ...  SEC 
NUMBER  OF  REGIONS  =  .  .  . 

RELATIVE  DIELECTRIC  CONSTANTS  =  .  .  . 

CONDUCTIVITIES  (MHO/M)  =  .  .  . 

SURFACE  BOUNDARIES  (CM)  =  .  .  . 

REGION  .  .  .  INTERIOR  POINT:  RADIUS  =  ...  CM  THETA  =  ...  DEG 
PHI  =  .  .  .DEG  ABSORBED-POWER  DENSITY  =  .  .  .  W/M**3 
REGION  .  .  .  INTERIOR  POINT:  RADIUS  =  ...  CM  THETA  =  ...  DEG 
PHI  =  .  .  .DEG  ABSORBED-POWER  DENSITY  =  .  .  .  W/M**3 


AVERAGE  ABSORBED-POWER  DENSITY  =  ...  W/M**3 
TOTAL  ABSORBED  POWER  =  .  .  .  WATT 
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noon 


PROGRAM  CSM 

ELECTROMAGNETIC  ENERGY  DEPOSITION  IN  A  CONCENTRIC 
SPHERICAL  MODEL  OF  THE  HUMAN  OR  ANIMAL  HEAD 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON  /COEFF/AN P, BNP, ALPNP, BETS P 

COMMON  FKP,BJNP,BHNP,CEX,BDP,P, DP, SIGP,EO, TIME, R, THETA, PHI 
INC, NORG, NMIN 

DIMENSION  BDP  (9)  ,SBDF(9)  ,EPSP(9)  ,SIGP(9)  ,  P  (1 0 1 )  , DP  (1 00) 
COMPLEX*1  6  FKP  (10)  ,  CFX ,  ANP  \1  000)  ,BNP(1000)  ,ALPNP(1000)  , 
1BETNP  (1000)  ,  B  JNP  (100)  ,  BHNP  (100)  ,  Z 
CALL  ERRSET  (208,0,-1,1) 

FIE=3. 14159265358 97 93D0 
RAD=1 80. DO/PIE 
EPSO=8.85416D-12 
VFL=2 .997924562D8 

C  ***  READ  CONTROL  PAF AMETEFS 

5  FEAD  (5,1C,END=110) FREQ, EO, TIME, STOPR , NORG , NOCF 
10  FOFMAT  (4E10. 0,215) 

C  ***  COMPUTF  COMPLEX  TIME  VARIATION 
OMEGA=2. D6*PIE*FREQ 
ARG=- OMEGA* TIME 

CEX=DCMPLX  (DCOS (AFG)  ,DSIN  (ARG) ) 

C  ***  READ  DIELECTPIC  PROPERTY  PARAMETERS 

READ  (5,20)  (EPSP(I)  ,SIGP(I)  ,I  =  1,NOFG) 

20  FORMAT  (8E10.0) 

C  ***  COMPOTE  COMPLEX  PFOPAG ATION  CONSTANTS 
FACl=OMSGA/VEL 
DO  30  1=1, NORG 
FAC2=EPSP (I) /2. DO 
FAC3=DSQRT (1.D0+ (1 .DO/(EPSO*OMEGA) **2)* (S IGP (I) /EPSP (I)  )  ** 
REKP=FAC1*DSQRT (FAC2*  (FAC3+1 . DO) ) 

FIMKP=FACl*DSQPT  (FAC2*  (FAC3-1 .DO) ) 

FKP  (I) =DCMPLX (REKP , FIMKP) 

30  CONTINUE 

FKP  (NORG  +  1) =DCMPLX (FAC1 , 0. DO) 

C  ***  READ  RADII  OF  SUFFACE  BOUNDARIES 

FEAD(5,20)  (SBDP(I) ,1  =  1, NORG) 

DO  35  1=1, NORG 
BDP  (I)  =SBDP  (I)/1.D2 
35  CONTINUE 

C  ***  PRINT  OUT  TITLE  AND  BASIC  INPUT  DATA 

WRITE  (6,40) FREQ, EO, TIME, NORG 

40  FOFMAT  (• 1  ELECTROMAGNETIC  ENERGY  DEPOSITION  IN  A  CONCENTRIC 

1AL  MODEL  OF  THE  HUMAN  OR  ANIMAL  HEAD* /« -FREQUENCY  =',F9.2, 
2  FIELD  STRENGTH  =',F7.2,'  V/M  TIME  =',F7.2,'  SEC 

3P  OF  REGIONS  =',I3) 

WRITE  (6,41)  (EPSP  (I)  ,1=1, NORG) 

41  FOFMAT  (' ORELATIVE  DIELECTRIC  CONSTANTS  = • , 9  (F7. 2 , 2X) ) 

WRITE  (6,42)  (SIGP(I)  ,1=1, NORG) 

42  FOFMAT  (• OCONDUCTIVITIES  (MHO/M)  =• , 9 (F7. 3 , 2X) ) 

WRITE  (6,43)  (S BDP  (I)  ,I=1,NOPG) 

43  FOPMAT  ('OSURFACE  BOUNDARIES  (CM)  =• , 9 (F7. 3 , 2X) ) 

C  ***  COMPUTE  SERIES  EXPANSION  COEFFICIENTS  FOR  ELECTRIC 

C  ***  FIELDS 

CALL  COEF (ANP, BNP, ALPNP,BETNP, NMIN) 

WRITE  (6,45) 

45  FOFMAT  (' O' ) 

DO  70  1=1 , NOCR 

C  ***  READ  DEFINING  CHARACTERISTICS  OF  INTERIOR  POINTS  AT 
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CSM0001 
CSM0002 
CSM0003 
CSM0004 
CSM0005 
CSM0006 
, STOPR,  CSM0007 
CSM0008 
CSM0009 
CSM0010 
CSM0011 
CSM0012 
CSM0013 
CSM0014 
CSM0015 
CSM0016 
CSM0017 
CSM0018 
CSM0019 
CSM0020 
CSM0021 
CSM0022 
CSM0023 
CSM0024 
CSM0025 
CSM0026 
CSM0027 
CSMC028 
CSM0029 
CSM0030 
2)  CSM0031 

CSM0032 
CSM0033 
CSM0034 
CSM0035 
CSM0036 
CSM0037 
CSM0038 
CSM0039 
CSM0040 
CSM0041 
CSM0042 
CSM0043 
SPHEBICCSM0044 
MHZ  CSM0045 
NUMBECSM0046 
CSM0047 
CSM0048 
CSM0049 
CSM0050 
CSM0051 
CSM0052 
CSM0053 
CSM0054 
CSM0055 
CSM0056 
CSM0057 
CSM0058 
CSM0059 
CSM0060 


n  o 


C  ***  WHICH  ABSORBED-POWER  DENSITIES  APE  TO  BE  COMPUTED 

HEAD  (5,50) NPEG,F,THETAD,PHID 
50  FORMAT  (I5,3E10.3) 

S AVR=P 
P=F/1.D2 

THETA=THETAD/P AD 
PHI=PHID/PAD 
Z=FKP  (NPEG) *B 

CALL  BJYH  (B JNP , BHNP , Z , NC , STOPE) 

NC=NCr  2 

IF  (NC.GT.NMIN) NC=NMI N 
CALL  PL (THETA, NC,P, DP) 

C  ***  ABSOEBED-POWEE  DENSITY  AT  GIVEN  POINT  INTEFIOB  TO 


CSM0061 
CSM0062 
CSM0063 
CSM0064 
CSM 0065 
CSM0066 
CSM0067 
CSM0068 
CSM0069 
CSM0070 
CSM0071 
CSM0072 


C  ***  ABSOEBED-POWEE  DENSITY  AT  GIVEN  POINT  INTEFIOB  TO  P-TH  EEGIONCSM0073 

CALL  EVEC  (NEEG, PD)  CSM0074 

PD=.5D0*SIGP  (NFEG) *PD  CSM0075 

C  ***  PRINT  OUT  PARTICULARS  OF  INTEFIOB  POINT  OF  REGION  P  CSM0076 

WRITE  (6,6  0)  NREG,SAVF,THETAD,PHID,PD  CSM  007"? 

60  FORMAT  ('  REGION', 12,'  INTERIOR  POINT:  RADIUS  =',F8.3,'  CM  THETA  =CSM0078 
1 ' , F7. 2, '  DEG  PHI  =',F7.2,'  DEG  ABSORBED  POWEF  DENSITY  = • , FI  2. 8, • CSM 0079 


C  *** 


C  *** 


2  W/M**3' ) 

70  CONTINUE 

NN=NOEG*NMIN 

FAC=2.D0*PIE/(FAC1*FAC1) 

QS=0. DO 
QT=0 . DO 

DO  90  N=1 , NMTN 
FACN=2. D0*N+1 . DO 

QT=QT+FACN*DEEAL (ALPNP (NN+N) +BETNP (NN+N) ) 

QS=QS+FACN*  (CDABS  (ALPNP (NN+N) ) **2  +  CDABS (BETNP (NN+N) ) **2) 
90  CONTINUE 

QA=FAC* (DABS (QT) -QS) 

**  TOTAL  ABSORBED  POWEF 

TOTPO  W=2 . 6544 ID- 3*EO**2*QA/2 .DO 
**  AVEPAGE  A B SO F BED- POWER  DENSITY 

P AVG=TOTPOW/  (4.DO*PIE*BDP(NORG)  **3/3. DO) 


CSM  0080 
CSM0081 
CSM0082 
CSM0083 
CSM  0084 
CSM0085 
CSM  0086 
CSM  0087 
CSM  0088 
CSM  0089 
CSM  0090 
CSM  0091 
CSM  0092 
CSM0093 
CSM  0094 
CSM0095 


C  *** 

c  *** 


PRINT  OUT  AVEFAGE  ABSORBED-POWER  DENSITY  AND  TOTAL  ABSORBED  CSM0096 


***  POWER  CSM0097 

WRITE  (6,100) PAVG,TOTPOW  CSM0099 

100  FORMAT ('O', 9X, 'AVEFAGE  A BSOF BED-POW ER  DENSITY  =',1PD13.5,'  W/M**3 • CSM0099 
1/*0’ ,9X, 'TOTAL  ABSORBED  POWER  =',D13.5,'  WATT')  CSM0100 

GO  TO  5  CSM 01 01 

110  STOP  CSM  01 02 

END  CSM  01 03 

SUBROUTINE  COEF (ANP , BNP, ALPNP , BETNP, NMI N)  CSM0104 

GENERATES  EXPANSION  COEFFICIENTS  CSM0105 

IMPLICIT  PEAL*8  (A-H,0-Z)  CSM0106 

COMMON  FKP, B JNP, BHNP, CEX,BDP,P, DP, SIGP,EO, TIME, F, THETA, PHI, STOPE,  CSM0107 
INC, NOBG  CSM0108 

DIMENSION  NTEF  (10)  ,BDP(9)  ,SIGP  (9)  ,P  (101)  ,DP  (100)  CSM0109 

COMPLEX*!  6  ANP  (1000)  ,  BNP  (1000)  ,  ALPNP  (1000)  ,  BETNP  (1000)  ,BJHP1  (1000)  CSM 01 10 
1 ,  BJHP2  (1000)  ,  BJNP  (100)  ,  BHNP  (100)  ,SJNP1  (100)  ,DELNP,SNT11 ,  CSM 01 11 

2SNT12 ,SNT21 ,SNT22 ,TNT1 1 ,TNT1 2 , TNT21 ,TNT22 , ETA PI , ETAP2 , ZEP1 , ZEP2 ,  CSM0112 

3SNP11 ,SNP12 ,SNP21,SNP22,TNP11,TNP12,TNP21,TNP22,DEL1,DEL2,FKP(10)  , CSM 01 13 
4CEX , F ATIO , S HNP1 (100) ,Z  CSM0114 

COMPUTE  EXPANSION  COEFFICIENTS  AN1 , BN1 , A NN, BNN , ALPN1  ,  BETN1 ,  CSM0115 

ALPNN , BETNN  CSM0116 

N 1  =  0  CSM  0117 

N2=0  CSM0118 

DO  15  NP=1 , NOFG  CSM0119 

Z=FKP  (NR) *BDP  (NR)  CSM0120 
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CALI  BJYH  (BJNP,BHNP,Z,N, STOPS) 

1 

CSM0121 

* 

DO  5  1=1, N 

CSM0122 

SJNP1  (I) =BJNP  (I) 

CSM0123 

5  SHNP1  (I)  =BHNP  (I) 

CSM0124 

Z=FKP  (NH+1) *BDP (NE) 

CSM0125 

CALL  BJYH  (B JNP , BHKP , Z , NN, STOPS) 

CSM0126 

NMTN=  MINO  (N , NN) 

CSM0127 

NTER (NF) =  NMIN 

CSM0128 

N2=N2+NMIN 

CSM0129 

DO  10  1=1  ,NMIN 

CSM01 30 

BJHP1  (N1 +1) =S JNP1  (I) 

CSM0131 

BJHP1  (N2  +  I) =SHNP1  (I) 

CSM0132 

BJHP2  (N1+I)  =BJNP  ( I) 

CSM0133 

BJHF2  (N2+I) =  BHNP  (I) 

CSB0134 

10  CONTINUE 

CSM0135 

N 1=N 1 +2*NMIN 

CSM0136 

N2=N2+NMIN 

CSM0137 

j 

15  CONTINUE 

CSM0138 

NMIN=MINN (NTEP , NOEG) 

CSM0139 

NMIN=NKIN-2 

CSM0140 

DO  17  1=1 ,  NMIN 

CSM0141 

ALPNP  (I) =DCMPLX (0.D0,0.  DO) 

CSM0142 

17  BETNP  (I) =DCMPLX (0- DO, 0. DO) 

CSM0143 

- 

NSUM=NOEG*NMIN 

CSM0144  j 

DO  30  1=1, NMIN 

CSM0145 

% 

JJ  =  0 

CSM0146  j 

KK  =  0 

CSM0147 

111=1+1 

CSM0148  j 

»> ' 

112=2*1+1 

CSM0149  1 

SNT 1 1 =DCMPLX (1 . DO , 0. DO) 

CSM015C  \ 

SNT1 2  =  DCMPLX  (0.D0,0.D0) 

CSM0151 

SNT21=SNT12 

CSM0152 

SNT22=SNT1 1 

CSM0153 

TNT1 1 =SNT1 1 

CSM0154 

TNT12=SNT12 

CSM0155 

TNT21=SNT12 

CSM0156  ! 

TNT22=SNT1 1 

CSM0157 

DO  27  J=1 ,NOFG 

CSM0158 

KK  =  KK  +  NTEE (J) 

CSM01  59 

i. 

ETAP1=  (III *BJHP1  (JJ  +  I) -I*BJHP1  (JJ  +  I  +  2) )/II2 

CSM0160 

ETAP2=  (II1+BJHP2 (JJ  +  I) -I*BJHP2 (JJ  +  I  +  2)  )  /II2 

CSM0161 

9 

ZEP 1=  (III *BJHP1  (KK+I) -I*B JHP1 (KK  +  I  +  2) ) /II2 

CSM0162 

ZEP2=  (II1*BJHP2 (KK+I) -I*BJHP2 (KK+I  +  2) ) /II2 

CSM0163 

DELNP  =  BJHP1  (JJ+I  +  1) *ZEP1-BJHP1 (KK  +  I  +  1)  *ETAP1 

CSM0164  ; 

E ATTO=FKP  (J  +  1 ) /FKP { J) 

CSM0165 

SNP1 1=  (ZEP1 *BJHP2  (JJ+I+1) -FATI0*BJHP1 (KK+I  +  1)  *ETAP2) /DELNP 

CSM0166 

SNP1 2=  (ZEP1 *BJHP2  (KK  +  I  +  1 ) -EATI0*BJHP1  (KK+I+1) *ZEP2) /DELNP 

CSM0167 

SNP21 =  (FATI0*BJHP1 (JJ  +  I  +  1) *ET AP2- ET API* BJHP2 (JJ  +  I  +  1 )  ) /DELNP 

CSM0168  j 

X 

SNP22=  (EATI0*BJHP1  (JJ  +  I  +  1) *ZEP2~ETAP1*BJHP2 (KK  +  I+1) ) /DELNP 

CSM0169 

Z=S  NT 1 1 

CSM0170  j 

« 

SNT11=SNT11*SNP11+SNT12*SNP21 

CSM0171 

SNT12=Z*SNP12+SNT12*SNP22 

CSM0172 

Z=SNT21 

CSM0173 

SNT21=SNT21*SNP11+SNT22*SNP21 

CSM0174 

SNT22=Z*SNP12+SNT22*SNP22 

CSM0175 

TNPl 1 =  (PATI0*ZEP1*BJHP2 (JJ  +  I  +  1) -BJHP1 (KK+I  +  1)  *ETAP2)  /DELNP 

CSM0176 

?NP12=  (PATI0*ZEP1*BJHP2 (KK+I  +  1) -BJHP1 (KK+I  +  1)  *ZEP2)  /DELNP 

CSM0177 

TNP21 =  (BJHP1  (JJ  +  I  +  1) *ETAP2-E ATI0*ETAP1*BJHP2 (JJ+I  +  1) ) /DELNP 

CSM0178 

TNP22=  (BJHP1  (JJ  +  I  +  1) *ZEP2-RATI0*ETAP1*BJHP2  (KK  +  I  +  1) ) /DELNP 

CSM0179 

Z=TNT1 1 

CSM0180 
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TNT11=TNT11*TNP11+TNT12*TNP21 

TNT12=Z*TNP12+TNT12*TNP22 

Z=TNT21 

TNT21=TNT21*TNP11+TNT22*TNP21 
TNT  22  =  Z  *TNP12+TNT22*TNP22 
JJ=JJ+2*NTEB  (J) 

KK=KK  +  NTEF  (J) 

27  CONTINUE 

ANP (I) =SNT1 1  - (SNT12*SNT21) /SNT22 
BNP (I) =TNT1 1- (TNT1 2*TNT21 ) /TNT22 
LL=NSUM+I 

ANP  (LL) =DCMPLX (1 .DO, 0. DO) 

BNP  (LL) =DCMPLX (1 .DO, 0.  DO) 

ALPNP  (LL) =-SNT21/SNT22 
BETNP  (LL) =-TNT21/TNT22 
30  CONTINUE 

IF  (NOKG.EQ. 1)  BETUFN 

C  COMPUTE  EXPANSION  COEFFICIENTS  AN 2 , . . . , AN  (N-1 )  ; BN2 , . . . , 

C  BN  (N-1)  ; ALPN2, . . . ,ALPN (N-1) ; BETN2, . . . ,BETN  (N-1) 

JJ=0 

KK=0 

MMl=0 

MM2=NMIN 

NFGM1=NOBG-1 

DO  45  J=1 , NBGM1 

KK=KK+NTEB (J) 

DO  40  1=1 , NMIN 
111=1+1 
112=2*1+1 

ETAP1=  (II1*BJHP1  (JJ  +  I)  -I*BJHP1  (JJ+I  +  2)  )  /II2 
ETAP2=  (II1*BJHF2 (JJ  +  I) -I*BJHP2 (JJ  +  I  +  2)  )  /II2 
ZEP1= (III *B JHP1  (KK  +  I) -I*BJHP1  (KK  +  I+2) ) /II2 
ZEP2=  (III *BJHP2 (KK  +  I) -I*BJHP2 (KK  +  I  +  2) ) /II2 
DELNP=BJHP1 (JJ+I+1) *ZEP1-BJHP1 (KK+I+1 ) *ETAP1 
FATIO=FKP (J+1) /FKP  (J) 

SNP1 1= (ZEP1*BJHP2  (JJ  +  I  +  1) -PATIO* BJHP1  (KK+I  +  1) *ETAP2) /DELNP 
SNP1 2= (ZEP1 *BJHP2 (KK  +  I  +  1 ) -EATIO*BJH PI  (KK+I  +  1) *ZEP2) /DELNP 
SNP21 = (F  ATI0*BJHP1 (JJ  +  I  +  1) *ETAP2-ETAPl *BJHP2 (JJ+I  +  1) ) /DELNP 
SNP22= (BATI0*BJHP1 (JJ+I+1) *ZEP2-ETAP1 *BJHP2 (KK+I+1) ) /DFLNP 
DEL1=SNP11*SNP22-SNP12*SNP21 

TNP11= (BATI0*ZEP1*BJHP2 (JJ+I+1) -BJHP1 (KK+I+1) *ETAP2) /DELNP 

TNPl 2=  (PATIO*ZEPl*BJHP2  (KK  +  I  +  1) -BJH Pi (KK+I  +  1) *ZEP2) /DELNP 

TNP21=  (BJHP1  (JJ  +  I  +  1)  *ETAP2-F.  ATI0*ETAP1*BJHP2  (JJ+I  +  1)  )  /DELNP 

TNP22=  (BJHP1  (JJ  +  I  +  1) *ZEP2-KATIO*ETAP1*BJHF2 (KK+I  +  1) ) /DELNP 

DEL2=TNP11*TNP22-TNP12*TNP21 

NN1=MM1+I 

NN2=MM2+I 

ANP (NN2)  =  (ANP (NN1) *SNP22- ALPNP (NN1 ) *SNP12) /DELI 
BNP (NN2)  =  (BNP (NN 1 ) *TNP22-BE?NP (NN1)*TNP12) /DEL2 
ALPNP  (NN  2)  =  (-ANP (NN1 ) *SNP21  + ALPNP  (NN 1 ) * SNP1 1 ) /DEL 1 
BETNP (NN2)  =  (-BNP (NN1 ) *TNP21  + BETNP  (NN1 ) *TNP1 1 ) /DEL2 
40  CONTINUE 

JJ=JJ  +  2*NTEE  (J) 

KK  =  KK  +  NTE F  (J) 

MMl=MM1 +NMIN 
MK2=MM2+NMIN 
45  CONTINUE 
FETUFN 
END 

SUBPOUTINE  EVEC  (NP, PD) 
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AO-AOa5  062  SCHOOL  Of  AEROSPACE  MEDICINE  BROOKS  AFB  TX  -  F/6  6/ia 

ELECTROMAGNET I C  ENERGY  DEPOSITION  IN  A  CONCENTRIC  SPHERICAL  MOD— ETC(U) 
DEC  79  EL  BELL,  D  K  COHOON,  J  M  PENN  TC 

UNCLASSIFIED  SAM-TR-79-6  w 


non 


COMPUTES  THE  FADIAL, COLATITUDE, AND  AZIMUTHAL  CSM0241 

COMPONENTS  OF  ELECTRIC  FIELD  VECTOE  E  FOP  CSM0242 

FEGION  P  AND  SCALAE  PRODUCT  E.E*  CSM0243 

IMPLICIT  EEA1*8  (A-H,0-Z)  CSM0244 

COMMON  /COEFF/ANP,BNP,ALPMP, BETNP  CSM0245 

COMMON  FRP,BJNP,BHNP,CEX,BDP,P,DP,SIGP,EO,TIME,B,THETA,PHI,STOPR,  CSM0246 
1NC,N0BG,NMIN  CSN0247 

DIMENSION  BDP(9)  ,SIGP(9)  ,P(101)  ,DP(100)  CSM0248 

COMPLEX+1 6  EF  AD, ETHETA , EPHI  ,FKP  (1  0)  ,ANP(1000)  ,  BNP  (1000)  ,  CSM0249 

1 ALPNP  (1000)  ,  BETNP  (1000)  ,BJNP  (100)  , BHNP  (100)  , CEX,T,Tl ,  PF.OD  CSM0250 

EBAD=DCMP1X  (0.D0,0.D0)  CSM0251 

ETHETA=DCMPLX  (0.D0, 0.D0)  CSM0252 

EPHI=DCMPLX  (0.D0/0.D0)  CSM0253 

NCK=0  CSM0254 

NN= (NP-1) *NMIN  CSM0255 

DO  40  N=1 ,NC  CSM0256 

FAC1=2. D0*N+1 . DO  CSM0257 

FAC2=N*(N+1.D0)  CSM0258 

E ATIO=FACl/FAC2  CSM0259 

T=FACl *P (N) *  (BNP  (NN  +  N) *BJNP  (N  +  1 ) + BETNP  (NN+N) *BHNP  (N  +  1) )  CSM0260 

NCR=NCR+1  CSM0261 

CALL  TEEM  (NCR,T,1)  CSM0262 

EFAD=EBAD+T  CSM0263 

T=ANP  (NN+N) *BJNP  (N  +  1 ) +ALPN? (NN+N) *BHNP (N+1)  CSM0264 

CALL  TERM  (NCR, T,0)  CSM0265 

T1=BNP (NN+N) * ( (N+1. DO) *BJNP (N) -N*BJNP (N+2) ) /FAC1+BETNP (NN+N) *  CSM0266 

1 ((N  +  1.D0)*BHNP(N)-N*BHNP(N*2))/FACl  CSM0267 

CALL  TEEM  (NCR, T1 ,1)  CSM0268 

IF ( (THETA. GE.1.D-6) .AND. (THETA.LT. 3. 141592D0) ) GO  TC  20  CSM0269 

IF(THETA.GE.3.141592D0)GO  TO  10  CSM0270 

ETHETA=ETHETA+FAC1/2.DO*T-RATIO*DP(N)*T1  CSM0271 

EPHI=EPHI-B ATIO+DP (N) *T+FACl /2 . D0*T1  CSM0272 

GO  TO  30  CSM0273 

10  ETHETA=ETHETA+(-1 .DO) ** (N+1) *FAC1/2 .DO*T-BATIO*DP (N) *T1  CSM0274 

EPHI=EPHI-EATIO*DP  (N) *T+  (-1 .DO) **  (N  +  1) *FAC1/2 . D0*T1  CSM0275 

GO  TO  30  CSM02'76 

20  ETHETA=ETHETA+FATIO*P (N) /DSTN (THETA) *T-FATIO*DP (N) *T1  CSM0277 

EPHI= EPHI- RATIO* DP  (N) *T+FATIO*P (N) /DSIN  (THETA) *T1  CSN0278 

30  IF(NCR.EQ.4) NCK=0  CSM0279 

40  CONTINUE  CSM0280 

PPOD=EO*CEX  CSM0281 

EFAD=-PHOD*DCOS (PHI) / (FKP (NP) *E) *EP AD  CSM0282 

ETHETA.  =  PPOD*DCOS  (PHI)  *ETHETA  CSM0283 

EPHI=PEOD*DSIN  (PHI) *EPHI  CSM0284 

PD=DPEAL  (EPAD*DCONJG (EFAD) ) +DBEAL  (ETHETA*DCONJG (ETHETA) ) +DBEAL (EPHCSMOI 


1 I+DCONJG (EPHI) ) 

CSN0286 

FETUPN 

CSM0287 

END 

CSM0288 

SUBROUTINE  TEEM (NCR, T,REY) 

CSM0289 

COMPUTFS  I**NCR*(N-TH  TEFM  IN  SERIES) 

CSM0290 

IMPLICIT  FEAL*8  (A-H,0-Z) 

CSM0291 

COMPLEX*1 6  T 

CSM0292 

IF (REY.EQ.I)GO  TO  20 

CSM0293 

GO  TO  (5,10,15,45) ,NCK 

CSM029U 

20  GO  TO  (10,15,45,5  )  ,  NCR 

CSM0295 

5  T  =  DCMPLX (-DIMAG(T)  ,DFEAL (T) ) 

CSM0296 

GO  TO  45 

CSM0297 

10  T=-T 

CSM0298 

GO  TO  45 

CSM0299 

15  T  =  DCMPLX (DIMAG (T)  ,-DFEAL  (T) ) 

CSM0300 
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45 

RETURN 

CSM0301 

END 

CSM0302 

SUBROUTINE  BJYH (BJNP, BHNP, Z , NN, STOPR) 

CSM0303 

c 

GENERATES  SPHERICAL  BESSEL  FUNCTIONS  JN  (KB)  AND  YN  (KR) 

CSM0304 

c 

AND  SPHERICAL  HANKEL  FUNCTIONS  OF  THE  PIPST  KIND  HN  (KR) 

CSM0305 

IMPLICIT  REAL*8  (A-H,0-Z) 

CSM0306 

COHPLEX*1 6  BJNP  (100)  , BYNP  (100)  ,BHNP  (100)  ,QP,Z 

CSM0307 

BYNP  (1)  =-CDCOS  (Z)/Z 

CSM0308 

BYNP  (2)  =  (BYNP  (1)  -CDSIN  (Z) )  /Z 

CSH0309 

DO  5  (1=3,100 

CSM0310 

BYNP  (M)  =  (2*M-3)/Z*BYNP(M-1) -BYNP  (M-2) 

CSH0311 

IF  (CDABS  (BYNP  (M) ) .GE. STOPR) GO  TO  10 

CSM0312 

5 

CONTINUE 

CSM0313 

10 

IF  ((!.  GT.  3)  GO  TO  25 

CSM0314 

c 

*** 

PRINT  OUT  ERROR  MESSAGE 

CSM0315 

WRITE  (6,20) Z 

CSM0316 

20 

FORMAT  (• 0****  PROCESS  CAN  NOT  PROCEED  SINCE  NM2  =  0 

FOR  Z  =' , 

1P2D1CSM0317 

15.7) 

CSM0318 

STOP 

CSH0319 

25 

BJNP  (M) =DCMPLX (0.D0,0.D0) 

CSM0320 

BJNP  (M-1)=-1.D0/(Z*Z*BYNP(M)  ) 

CSM0321 

NM2=M-2 

• 

CSM0322 

DO  30  1=1, NM2 

CSM0323 

L=M-I 

CSM0324 

BJNP  (L-1 )  =  (2*L-1 ) /Z*BJNP (L) -BJNP  (L+ 1) 

CSM0325 

30 

CONTINUE 

CSM0326 

QP=CDSIN  (Z) /  (Z*BJNP (1 ) ) 

CSM0327 

NM1=M-1 

CSM0328 

DO  35  N=1 , NMl 

CSM0329 

NN=N 

CSM0330 

BJNP  (N) =QP*BJNP (N) 

CSM0331 

IF  (CDABS  (BJNP  (N) ) . LT .1.D-25) GO  TO  40 

CSM0332 

35 

CONTINUE 

CSM0333 

40 

DO  45  1=1, NN 

CSM0334 

REJN=DFEAL  (BJNP (I)) 

CSM0335 

FIHJN  =  DIMAG  (BJNP (I) ) 

CSM0336 

FEYN=DBEAL (BYNP (I) ) 

CSM0337 

FTMYN=DIMAG  (BYNP (I) ) 

CSM0338 

PEHN=REJN-FIMYN 

CSM0339 

FIMHN=FEYN+FIMJN 

CSM0340 

BKNP  (I) =DCMPLX (REHN, FIHHN) 

CSM0341 

45 

CONTINUE 

CSM0342 

PETUPN 

CSM0343 

END 

CSMC344 

SUBPOUTINE  PL  (THETA, N,P, DP) 

CSM0345 

c 

GENERATES  ASSOCIATED  LEGENDRE  FUNCTIONS  OF  THE 

FIRST 

CSM0346 

c 

KIND,  OFDER  1  AND  DEGREE  N,  AND  THEIR  FIRST  DERIVATIVES 

CSM0347 

IMPLICIT  FEAL*8  (A-H,0-Z) 

CSM0348 

DIMENSION  P  (101)  ,DP(100) 

CSM0349 

SNJ=DSIN  (THETA) 

CSM0350 

CNJ=DCOS (THETA) 

CSM0351 

P  (1 )  =SNJ 

CSM0352 

P  (2)=3.D0*SNJ*CNJ 

CSM0353 

DP  (1)  =CN  J 

CSM0354 

DO  10  M=2,N 

CSM0355 

A=M 

CSM0356 

MP1=M*1 

CSM0357 

P  (MP1)  =(2.D0*A+1.D0)/A*CNJ*P (M) - (A+1.D0)/A*P (M-1) 

CSM0358 

IF  ( (THETA. GE.1.D-6) .AND. (THETA . LT.  3. 1 4 1 592D0) )  GO  TO 

5 

CSM0359 

DP  (M) =M*MP1/2 

CSM0360 
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r 


IF  (THETA .GE.3.141592D0)  DP  (M)  =  (-1  .DO)  **M*DP  (M) 

CSM0361 

GO  TO  10 

CSM0362 

5 

DP(M) =  (A*CNJ*P  («)- (A+1.D0)  *P (M-1) )/SNJ 

CSH0363 

10 

CONTINUE 

CSM0364 

RETUPN 

CSM0365 

END 

CSM0366 

FUNCTION  MINN  (NB AY, N) 

CSK0367 

DETERMINES  MINIMUM  POSITIVE  INTEGEP  VALUE 

CSM0368 

DIMENSION  NPAY{10) 

CSM0369 

IF  (N. EQ. 1 ) GO  TO  20 

CSM0370 

NMIN=NRAY (1) 

CSM0371 

DO  10  1=2, N 

CSM0372 

NTEMP=NP AY (I) 

CSM0373 

IF  (NTEMP. LT. NHIN) NMIN=NTEMP 

CSM0374 

10 

CONTINUE 

CSM0375 

MINN=NMIN 

CSM0376 

GO  TO  30 

CSH03",7 

20 

H INN  =  NF A  Y  (1) 

CSM0378 

30 

EETUPN 

CSB0379 

END 

CSM0380 
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